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Preface 



The main question that this thesis tries to answer is the following: How can we calculate 
numerically the gravitational field of isolated systems in an efficient way based upon a 
clear geometric framework? 

Isolated systems in general relativity are models of self-gravitating astrophysical 
sources. Depending on the scales in consideration, these sources can be planets, stars, 
black holes and even galaxy clusters. It is not the source that makes a system isolated, 
but the behavior of its gravitational field far away from the source. The systems that 
we might consider isolated share the property that as we move away from the source, 
the gravitational field becomes weak. We model such systems by suitably attaching to 
the far-field zone an asymptotic region in which the spacetime becomes flat in a certain 
sense. These models are commonly referred to as asymptotically flat spacetimes. 

In this idealization procedure, the details of the asymptotic behavior need to be 
decided upon in accordance with physical intuition and the equations describing the 
model, as we can not yet observe the precise behavior of asymptotic gravitational fields 
to test and modify our assmnptions. 

An elegant and powerful characterization of the asymptotic behavior of spacetimes 
can be given in terms of the conformal structure. Conformal techniques make use of 
an important interplay between physics, conformal geometry and the theory of partial 
differential equations in the study of Einstein equations. The conformal approach has 
been applied successfully in mathematical relativity to various problems concerned 
with the large scale structure of spacetimes. 

To relate observations of gravitational radiation to properties of astrophysical 
sources, we need to understand spacetimes representing isolated systems arising from a 
large class of initial data that correspond to astrophysical configurations. As the Ein- 
stein equations can not be solved explicitly for interesting dynamical systems emitting 
gravitational radiation, we need to apply approximation schemes to make quantita- 
tive predictions that can be compared with observations. The post-Newtonian and 
perturbative schemes cover the weak field and low velocity limits. To make reliable 
predictions on highly dynamical strong fields, we need numerical simulations. 

There are many open problems in numerical calculations of spacetimes. At least 
two of these problems, namely the treatment of the outer grid boundary and the inter- 
pretation of numerically generated spacetimes, are related to the asymptotic region. 
Applying conformal techniques to deal with these problems can lead to accurate and 
efficient codes that allow us to gain control over numerical approximations. 

An important goal of this thesis is the development of a numerical code that would 
allow us to calculate entire, asymptotically flat, radiative solutions to the Einstein 
equations. With such a code we should be able to follow the maximal development of 
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asymptotically flat initial data starting from a Cauchy surface reaching spatial infinity 
up to the region close to timelike infinity. The achievement of this goal seems to 
require different methods adapted to different asymptotic regions, i.e. spatial infinity, 
null infinity and timelike infinity. With this goal in mind, we will discuss two novel 
applications of conformal techniques to numerical simulations. These new methods 
share the important advantage that the representation of the conformal factor is known 
a priori in terms of grid coordinates. I will argue that this feature is very convenient, 
because it is crucial to have numerical access to and some control over the asymptotic 
region for accurate and efficient calculations of asymptotically flat spacetimes. 

Organization 

In the introduction we describe shortly the development of some basic notions in con- 
formal geometry. We mention the controversy on the nature of gravitational radiation 
and then present the idea of conformal infinity introduced by Penrose that is widely 
used in mathematical relativity today. We discuss conditions for the feasibility of the 
conformal approach investigated by Friedrich. An overview of current methods in nu- 
merical relativity with respect to the treatment of the far-field region of asymptotically 
flat spacetimes is given to motivate the numerical application of conformal techniques 
further. 

Chapter [2] concentrates on infinity in null directions, called null infinity. In l2.1l and 
12.21 we investigate spacelike slices reaching null infinity, so-called hyperboloidal slices, 
and conformal compactification of the Minkowski and the extended Schwarzschild 
spacetimes. We see that one can prescribe the representation of a conformal factor in 
terms of a suitably chosen compactifying radial coordinate on spherically symmetric 
hyperboloidal slices such that null infinity is at a fixed spatial coordinate location. 

In l2.3l and l2.4l we construct a method to numerically treat the hyperboloidal initial 
value problem including null infinity in the computational domain without making 
any symmetry assumptions. This method allows us to find solutions to the Einstein 
equations in coordinates in which null infinity is fixed to a spatial coordinate location. 
It is based on the general wave gauge that is commonly used in numerical relativity. It 
introduces a suitable coupling of the conformal and the coordinate gauge to guarantee 
certain geometric properties of null infinity with an appropriate choice of the gauge 
source functions for the coordinates. It requires the numerical calculation of formally 
singular terms arising from the conformal compactification. 

In 12.51 we discuss numerical test results obtained with the suggested method in 
the special case of spherical symmetry on the example of the extended Schwarzschild 
spacetime. We use a simple choice of evolution variables and numerical boundary 
treatment. While our numerical setup does not allow us to do long time evolutions 
of the extended Schwarzschild spacetime, one can see that the method can be applied 
even with crude numerical techniques without the formally singular terms leading to 
an immediate blow-up of numerical errors at the outer boundary. 

The chapter ends in l2.6l with an outlook on next steps to establish this new approach 
for simulating dynamical isolated systems, and a short discussion of its open problems 
and limitations. 

In chapter[3]we include not just null infinity but also spatial infinity in the numerical 
domain which gives us direct access to the global structure of asymptotically flat 
spacetimes. The reduced general conformal field equations developed by Friedrich 
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provide the only available system for the treatment of spatial infinity. The system is 
based on the conformal Gauss gauge that we study using numerical methods. In l3.f I we 
reproduce Friedrich's construction of a conformal Gauss gauge in the Schwarzschild- 
Kruskal spacetime numerically covering the entire solution in a smooth way. Going 
beyond analytical studies we find out numerically that one can also cover the Kerr 
solution using conformal geodesies including null infinity, timclike infinity and the 
Cauchy horizon. 

In 13.21 we discuss certain aspects of the reduced general conformal field equa- 
tions relevant for their numerical implementation. In 13.31 we solve numerically the 
Cauchy problem for the equations in spherical symmetry with initial data from the 
Schwarzschild-Kruskal spacetime which gives us the first numerical calculation of an 
entire asymptotically flat black hole spacetime including timelike, null and spacelike 
inflnity and the region close to the singularity. 

In order to include dynamical gravitational fields into our discussion of spacetimes 
in a neighborhood of spatial infinity, we calculate in 13.41 asymptotically fiat, axially 
symmetric initial data based on studies by Friedrich. The calculated initial data has the 
special property that its ADM-mass vanishes but its development has a non- vanishing 
radiation field. In l3.5l we develop this data on a three dimensional Cartesian grid using 
the reduced general conformal field equations such that spatial infinity is represented 
by the point at the origin of our Cartesian grid. We calculate a certain component 
of the Weyl tensor in a suitably adapted Newman-Penrose tetrad representing the 
radiation field at null infinity and show that it does not vanish along null infinity in 
accordance with expectation. 

The Cartesian code can not be used to study spacetimes with non- vanishing ADM- 
mass. For a generalizable code we implement the regular finite initial value problem 
near spatial infinity formulated by Friedrich based on the reduced general conformal 
field equations in a gauge that allows us to represent spatial infinity as a cylinder. 
The cylinder at infinity imposes a certain geometry that requires a numerical code 
which can handle a spherical grid topology in a frame formalism. In 13.61 we discuss 
a numerical implementation of a frame-based evolution system using spherical grid 
topology in three spatial dimensions with overlapping grids. 

In 13.71 we implement the regular finite initial value problem with the cylinder at 
spatial infinity. A difficulty caused by the degeneracy of the equations at the set 
where the cylinder at spatial infinity touches null infinity is dealt with by freezing 
the evolution in the unphysical domain by exploiting the a priori knowledge of the 
conformal factor in terms of initial data and grid coordinates. The chapter ends with 
a discussion of possible improvements of the code. 

In chapter [4] I summarize the main results of the thesis and give an outlook for 
possible directions for future work. The chapter ends with remarks on the idea of 
conformal infinity. 

In the appendix we discuss the generation of causal diagrams. 



Chapter 1 

Introduction 



This thesis is concerned with numerical studies of asymptotic fields of isolated systems. 
Conformal techniques play a fundamental role in the discussion of the asymptotic struc- 
ture of spacetimes. To develop an intuitive understanding of conformal rescalings, I 
will refer to their application in cartography where their advantages and disadvantages 
can be addressed in a simple manner and describe the argumentation of Hermann Weyl 
to emphasize the role of the conformal structure of spacetimes in general relativity. 

The historical controversy concerning the physical nature of gravitational radia- 
tion is especially demonstrative for being cautious about interpretations of coordinate 
dependent calculations in general relativity. I will indicate how an unambiguous for- 
mulation of gravitational radiation that uses conformal techniques has been achieved. 

The Penrose conformal compactification technique will be discussed and results 
on its feasibility due to Friedrich will be described. It is a goal of this thesis to 
combine a rigorous treatment of the far-field region of asymptotically flat spacetimes 
with advantages of widely used numerical methods. To motivate this goal I will point 
out certain properties of common numerical methods with special emphasis on their 
treatment of the asymptotic region. 

The introduction follows chronological lines because I believe that some of the 
current problems in numerical relativity have their counterparts in earlier problems 
that have been thoroughly discussed. We can learn from these discussions to avoid 
similar confusions. Naturally, the presentation is highly subjective and makes no claim 
of completeness. 

1.1 Conformal geometry 

A map that preserves local angles is called conformal. The first problem where con- 
formal mappings played an important role, beside the trivial rescaling where the con- 
formal factor is just a positive constant, was the problem of projecting the surface of 
a globe onto a plane. The stereographic projection, probably already known to the 
Egyptians, has been applied in the 2nd century BC by Hipparchus to represent the 
celestial sphere on an astrolabe. Ptolemy used the stereographic projection for his car- 
tography, a field which historically saw the most influential application of conformal 
techniques [57]. The Mercator projection of the Earth constructed in the 16th century 
is also a conformal projection. It was devised for nautical charts with the property 
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that all lines of constant compass direction are represented by straight segments, but 
its use was not only restricted to marine navigation. The Mercator projection was 
so influential that its wide use was subject to political controversies related to the 
distortion of scales in conformal maps [38j . 

The development of methods in cartography has paved the way for the discussion 
on geometry of curved surfaces by Lambert, Bolyai, Lobachevsky, Gauss, Riemann and 
others. Indeed, some of the basic terminology in differential geometry has its origin 
in cartography. Therefore, it does not surprise that certain problems of cartography 
are similar to certain problems of general relativity. For example, it is clear today 
that in cartography different map projections serve different purposes. There is no 
single map of the Earth that covers all requirements. Similarly, there is no best gauge 
choice for a spacetime in general relativity. Emphasis on different aspects and parts 
of the solution require different gauge choices as will be demonstrated in this thesis 
by studies of various representations of the Minkowski and the Schwarzschild-Kruskal 
spacetimes. 

A key feature leading to the wide use of the Mercator projection is the fact that it is 
a conformal projection. As a consequence the projection keeps relative local directions 
-in the sense of local angles- invariant. Clearly, relative local directions are important 
for marine navigation. In general relativity, the fundamental important structure is 
the causal structure which is equivalent to the null cone structure as well as to the 
conformal structure. We will see in the next sections how this fact can be used to get 
access to the large scale structure of spacetimes and to describe gravitational radiation 
rigorously, in an unambiguous way. 

The importance of conformal geometry in general relativity has been recognized 
already in 1918 by Hermann Weyl ,146,. His motivation to study conformal rescalings 
of the metric has partly been philosophical |126j . He was disturbed by what he called 
an inconsistency due to which lengths of vectors at different points in Riemannian 
geometry can be directly compared with each other. Weyl observes in |146j : "The 
metric allows the magnitudes of two vectors to be compared, not only at the same 
point, but at any two arbitrarily separated points." To remove this remnant of non- 
local geometry (" Ferngeometrie" ) , he proposed to consider the conformal class of a 
spacetirne^ 

Let M he a, four dimensional, smooth manifold and g a Lorentzian metric. The 
conformal class [g] of a spacetime {A4,g) is given by metrics g related to g by a 
conformal rescaling g — fl^g with a positive, point dependent conformal factor > 0. 
A choice of the unit of measurement, i.e. the gauge, corresponds to the choice of a metric 
in the conformal class [g] or equivalently to the choice of a point dependent conformal 
factor n. As a consequence, in Weyl's "pure infinitesimal geometry" |147j . the choice 
of the unit of measurement is subject to local variations. Weyl states in [146j that for 
the explicit representation of the spacetime in such a geometry we have to choose a 
coordinate system and at each point determine the conformal factor. Therefore each 
formula must have a double-invariance, namely invariance with respect to arbitrary 
smooth coordinate transformations and with respect to conformal rescalings of the 
metric. For applications of conformal techniques in this thesis, especially in chapter O 
a suitable coupling of the coordinate and the conformal gauge freedom will be essential 
to gain control over the asymptotic region. 
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The main step from Euclidean geometry to Riemannian geometry is the removal 
of the assumption of integrability of vectors by parallel transport. In Riemannian 
geometry, a vector parallel transported along a closed curve changes in general its 
direction due to curvature, while its norm stays the same. Weyl demanded "the 
non- integrability of the transference of distances" [150j . He states in [146] : true 
infinitesimal geometry should, however, recognize only a principle for transferring the 
magnitude of a vector to an infinitesimally close point and then, on transfer to an 
arbitrarily distant point, the integrability of the magnitude of a vector is no more to 
be expected than the integrability of its direction." To allow for comparison of lengths 
in each small neighborhood, he introduced what we call today the Weyl connection. 
A concise description of these ideas can be found for example in [U I105j . We use a 
modern representation following the notation of 64J. 

A Weyl connection V is a torsion free connection which is not necessarily the Levi- 
Civita connection of a metric in the conformal class [g] but preserves the conformal 
structure so that the covariant derivative of any metric g in [g] is proportional to itself, 

Va5pi^ = -'^fxQtJ.v, (1-1) 

with a one-form /. Under conformal transformations of the metric g' — uj'^g with a 
function w > 0, the one-form / transforms according to /' = / — uj~^ duj. If / is 
exact, it can be written as f — il.~^dfl and the Weyl connection V is the Levi-Civita 
connection V of a metric g in the conformal class given hy g = fl^g. 

The relation of the Weyl connection V with a metric connection V of a metric g 
in the conformal class is given by 

V-V = 5(/), where S{f); , = S^J, + S^j'^ - ~g,,,r^h. (1.2) 

We see that Weyl connections characterized by 1-forms / allow us a more general 
discussion of the conformal structure than conformal rescalings of the metric. 

The introduction of the Weyl connection was the first step in the study of gener- 
alized connections and it led to a rich geometry, which Weyl used in his attempt to 
unify the gravitational and electromagnetic fields. His attempt for a unified theory 
was unsuccessful due to physical observations [76], however the main idea has been 
the basis during the construction of modern gauge theories |101] . It is interesting to 
note that the concept of gauge invariance in theoretical physics originates in Weyl's 
ideas on conformal invariance 102J. 

A torsion free connection implies the notion of parallel transport and allows the 
construction of geodesies. Auto-parallel curves with respect to a Weyl connection are 
referred to as conformal geodesies [TOj. Remarkably, conformal geodesies have been 
useful in investigations of the asymptotic behavior of solutions to Einstein equations 
carried out by Helmut Friedrich [64]. Following his investigations, we will make use of 
conformal geodesies in chapter [3] for a numerical treatment of spatial infinity. 

The application of conformal techniques in general relativity played an important 
role in the clarification of a heated controversy in mathematical relativity on the phys- 
ical nature of gravitational radiation that we discuss in the next section. 
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1.2 Gravitational radiation and asymptotic flatness 

The Newtonian theory of gravitation is governed by a hnear eUiptic Poisson equation 
for the Newtonian potential and does not allow for free gravitational degrees of freedom 
carrying energy. The action-at-a-distance property of Newtonian gravitation related 
to the elliptic nature of the underlying Poisson equation was already questionable to 
Newton and his contemporaries. A finite speed of gravity has been studied by Laplace 
who concluded from this assumption to a damping force on the orbital motion of 
planetary systems. In 1908, Poincare suggested an acceleration of planetary orbits 
due to loss of energy by emission of gravitational waves. 

The first relativistic description of gravitational waves is due to Albert Einstein [iSl 
147) . He studied wave phenomena in linearized gravity and introduced the quadrupole 
formula which, in the weak field limit, relates the motion of a source to generation of 
gravitational waves. In his linearized calculations he made use of a coordinate depen- 
dent quantity which he called the pseudo energy momentum tensor to calculate the 
energy flux of gravitational waves. Einstein's use of the coordinate dependent pseudo- 
tensor was criticized by many scientists of his time [31] , among others Eddington who 
suggested to put more emphasis on the curvature tensor, showed that certain classes 
of wave-like solutions to the linearized Einstein equations are not physical and pointed 
out limitations of linearized calculations |45j . 

Until the 1960s the status of gravitational waves was quite uncertain. Even Einstein 
seems to have believed for some time that gravitational waves did not exist [86] . The 
use of linearized calculations whose qualitative predictions are not necessarily valid 
in the full theory, the application of coordinate-dependent methods and the non-local 
nature of gravitational field energy have been sources of confusion. The difficulty with 
a local energy concept for the gravitational field lies in the background independence 
of general relativity. Due to the equivalence principle, the local gravitational field can 
be transformed away by passing to a freely falling frame of reference. 

With the work of Lichnerowicz on mathematical foundations of general relativity 
|94j it turned out that certain calculations supporting the view that gravitational waves 
are unphysical put too restrictive conditions on coordinates. Based on the development 
of mathematical foundations and a stronger emphasis on curvature invariants, the 
question on the status of gravitational waves was gradually answered by the works 
of Bondi, Newman, Penrose, Petrov, Pirani, Robinson, Sachs and Trautman among 
others. 

It is difficult to give a precise definition to radiation in general relativity as the 
gravitational field acts as its own source so that there is no background we can rely 
on. However, the gravitational field becomes weak far away from sources. When we 
concentrate our studies on an astrophysical system we assume that the interaction of 
this system with the rest of the universe to be negligible. In the study of such "isolated 
systems" the astrophysical system is modeled by an asymptotically fiat spacetime. This 
corresponds to an idealized situation in which massive sources are confined to spatially 
compact domains and the spacetime metric is required to approach a flat metric at 
infinity. One might hope to find a precise notion of radiation for such systems in the 
asymptotic region using the flat limit metric as a kind of background structure. 

For the discussion of gravitational radiation, it is natural to require asymptotic 
flatness in null directions. Historically, however, first discussions of asymptotic flatness 
have been made with respect to spatial infinity. Note that asymptotic fiatness is 
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an assumption that is " put in by hand" . Today there are several mathematically 
inequivalent definitions of the concept and we will refer to some of them in later 
sections. In essence, one imposes some condition determining the approach of the 
metric to a flat metric at (spatial or null) infinity. A formal definition based on 
the existence of distinguished classes of coordinates^^say near spatial infinity, can be 
roughly given as follows. We may call a spacetime {Ai,g) asymptotically flat at spatial 
infinity when there exists a coordinate system {i'^} in a neighborhood of spatial infinity 
with respect to which 

9fiv = fjpi.v + Ok{f~^), as r ^ oo, for some fc,p > 0, (1.3) 

where = Safix^x^ with a, (3 = 1, 2, 3. We denote by fj the Minkowski metric. Sap is 
the Kronecker symbol and the terms Ok{r~^) behave like 0(f~^"*"^) under differentia- 
tions of order j < k. The strength of the fall-off behavior is determined when we give 
a precise description of Ofc(r~^). A question is how strong this fall-off behavior needs 
to be so that suitable physical notions of total energy, linear and angular momentum 
can be defined. It is known, for example, that for p > 1/2 and fc > 1 a well-defined 
notion of total energy at spatial infinity can be associated with an asymptotically fiat 
spacetime |130) . 

The linearized theory suggests that gravitational perturbations travel along null 
directions. Therefore, to define gravitational radiation, one takes limits to infinity in 
null directions. One may define a notion of asymptotic flatness in null directions by 
requiring (jl.3p in a neighborhood of null infinity with respect to a suitable coordinate 
system. We might intuitively imagine that the spacetime should become simply close to 
the Minkowski spacetime in whatever direction we approach infinity. The geometric 
picture, however, is quite different and it is a delicate question how the notions of 
asymptotic flatness at spatial and null infinity fit together. We will come back to this 
issue at the end of the section. Let us flrst follow the historical development of the 
concept of gravitational radiation. 

An important step in the clarification of questions regarding gravitational radiation 
was the demonstration that outgoing gravitational radiation carries positive energy 
away from an isolated system [2TJ I120j . This analysis used a special class of coor- 
dinate systems which is natural for studying radiation phenomena. The coordinates 
are generated by distinguished families of outgoing null hypersurfaces. On the hy- 
persurfaces the luminosity distance f has been defined. Today such coordinates are 
called Bondi coordinates. The field equations have been studied along null directions 
and certain conditions at infinity have been put that are similar to the no incoming 
radiation conditions suggested by Trautman |138) . These conditions included the re- 
quirement that metric components have certain smoothness properties in 1/r in the 
limit f ^ oo along null directions. An important open problem was whether and to 
what extend these conditions excluded interesting solutions |120| . 

Instead of putting conditions on metric components corresponding to asymptotic 
flatness as above, it is more convenient to impose restrictions on certain components 
of the Weyl tensor C^^^^ representing the free gravitational fleld. Sachs analyzed the 
vacuum Bianchi identities VpC^^^ = and suggested that the Weyl tensor satisfies 
a "peeling property" along outgoing null geodesies relating the fall-off behavior of its 
components to the Petrov classification |119j. Shortly after this analysis, Newman and 
Penrose found a way to formulate fall-off conditions in an elegant way and showed that 
they imply the peeling behavior |100j . They introduced the spin frame formalism and 
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analyzed expansions of gravitational fields in 1 jf where f is an afline parameter along 
outgoing null geodesies. Imposing a certain uniform smoothness condition, they found 
out that if a particular complex tetrad component of the Weyl tensor with respect 
to a distinguished tetrad at infinity has an asymptotic behavior i/jq = 0(r^'^) and 
5^00 — 0{r^^) as f oo, one can derive the Sachs peeling behavior tpk — 0(f'"'~^) as 
f ^ oo for fc = 0, 1, 2, 3, 4. In analogy to electromagnetism, the field with the weakest 
fall-off in the order 1/f, namely "04, is interpreted as representing the outgoing radiation 
field. Note that this interpretation of -04 is strongly related to the Sachs peeling 
behavior and requires a distinguished choice of a tetrad field at infinity. In contrast 
to electromagnetism, where the linearity of the Maxwell's equations allows a clean 
global separation between near field terms and far field terms, in general relativity no 
such separation can be expected in general as the peeling behavior is valid only in the 
asymptotic region. In the strong field regime one would expect a mixture of the fall-off 
behavior. Therefore it is misleading to talk about 04 calculated with respect to some 
tetrad field in the interior as the radiation field. 

In the studies just discussed, the close relationship between gravitational radiation 
and asymptotic behavior of solutions to the Einstein equations has been demonstrated. 
However, the studies still relied on the use of special classes of coordinate systems and 
included awkward limits to infinity. Such limits arise because a rigorous meaning to 
gravitational radiation can only be given in the asymptotic region. 

To study physical properties of isolated systems, one would like to have coordinate 
independent definitions. The coordinate dependence in the early studies of gravita- 
tional radiation blurred the global geometric structure of asymptotically fiat space- 
times. In the early development of the field, the global picture was by no means clear, 
not even to those who contributed substantially to the development of the picture. 
Sachs, for example, refers to null infinity as "the cylinder at spatial infinity" in [120j . 

To motivate a coordinate independent approach to the concept of gravitational 
radiation further and to point out a subtlety regarding the notion of asymptotic fiatness 
in spacelike and null directions, we shall discuss an illustrative example taken from 
|110| . The Schwarzschild metric with mass m in standard Schwarzschild coordinates 
{t, fs,'&, reads 

~gs = -(l- dt^ ^i}^ ^ ^ ^"^'^^ 

where dcr^ = dd'^ -\- siv? d d^p^ . This metric obviously approaches the flat Minkowski 
metric f]s given by fjs = —dt^ + dfg + f^da^, as fg oo, so that we may call the 
Schwarzschild spacetime asymptotically flat at spatial infinity. This, however, does 
not by itself imply that the Schwarzschild spacetime is asymptotically flat in null 
directions. With regard to null geodesies the Schwarzschild metric (|1.4p differs greatly 
from the Minkowski metric ffs at large distances. Consider outgoing null geodesies 
given by u const., 1} = const., (p = const., where u is the Schwarzschild retarded 
time u~t—{fs + 2m ln(rs — 2m)). We see that as ^ 00 the value of the Minkowski 
retarded time t — fg is unbounded above along outgoing Schwarzschild null geodesies 
(the argument applies also for ingoing null geodesies) . Schwarzschild null geodesies do 
not correspond at all to Minkowski null geodesies with respect to fjs for large fs- 

The solution to this apparent problem is to relate a Minkowski metric to the 
Schwarzschild metric in a different way. Define a coordinate f* = -I- 2m ln(fs — 2m). 
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We can choose on the Schwarzschild manifold another Minkowski metric given by 
fj^ = —dt^ + dfl + f1 d(T^ . Now the described difficulty does not arise because outgoing 
null geodesies given hy u — t — f^ — const, have the same form for both metrics gs and 
77*. While the condition for asymptotic flatness (|1.3p has the same form in spacelike 
and null directions, its geometric interpretation, i.e. the way we relate the Minkowski 
metric fj to the curved spacetime metric g, is different. 

We mentioned that for studying gravitational radiation one might use the flat 
limit metric as a kind of background structure. This example makes clear that the 
relation of the flat limit metric to the curved spacetime is determined by the null cone 
structure. The null cone structure, however, is part of the unknown in the Einstein 
equations. It is therefore a difflcult question how the flat background structure can be 
used conveniently in dynamical spacetimes to discuss their asymptotic properties. 

It turns out that the structure that we are looking for can be given in terms of the 
idea of conformal infinity presented in the next section. This construction also allows 
to replace asymptotic calculations by local differential geometry. 

1.3 Conformal infinity 

Conformal techniques in general relativity played an important role in the investiga- 
tion of the solution space to the Einstein equations. The analysis of elliptic constraint 
equations implied by the Einstein equations on spacelike hypersurfaces has made in- 
tensive use of conformal methods [93l I153j . In the following, we will concentrate on 
conformal rescalings of the full spacetime metric. We will also ignore the application 
of conformal techniques in cosmological models (see [BDl for an overview and [T7] for 
a recent work). 

It was discovered by Penrose that a suitable notion of asymptotic flatness can be 
formulated in an elegant way using the conformal equivalence class of the metric. The 
conformal compactification technique enabled a geometric formulation of the fall-off 
behaviour. Penrose could also deduce the peeling behavior from a few assumptions on 
the conformal structure of the metric not favoring any coordinate systems |107[ 1108] . 
Since then, conformal techniques to study in detail the structure of gravitational fields 
in the asymptotic region in a coordinate independent way have been very fruitful in 
mathematical relativity. Many physical concepts like mass, momentum or gravitational 
radiation of isolated systems have been unambiguously defined in the asymptotic region 
of an isolated system using conformal techniques. 

Two important observations suggest a conformal approach for discussing radiation. 
The first observation is the importance of the null cone structure. Conformal rescalings 
of the metric preserve the null cone structure, which is equivalent to the characteristic 
and the causal structures. It turns out that essential features of free gravitational 
fields strongly related to the null cone structure can be discussed with respect to 
the conformal equivalence class of the physical metric. The second observation is 
related to the non-local feature of gravitational radiation which requires limits to 
infinity as discussed in the previous section. Conformal rescalings can be used to 
avoid such limits by compactifying the spacetime such that infinity corresponds to a 
finite hypersurface. Asymptotic calculations can then be replaced by local differential 
methods under certain conditions. 



1.3 Conformal infinity 



12 



The idea of conformal compactification is similar to the construction of the Riemann 
sphere in complex analysis by adjoining infinity to the complex plane. Due to the 
Lorentzian signature of the metric, infinity for a spacetime corresponds not to a point 
but has a richer structure, as Weyl already observed in |148j . To demonstrate the 
technique we discuss the conformal compactification of the Minkowski spacetime (see 
also [5Tl [Tn9l[T44| V 

The Minkowski metric fj in coordinates adapted to spherical symmetry is given by 

f] = ^dP + df'^ dcr^, on i e R, f > 0, (1.5) 

where da^ is the standard metric on the unit sphere S^. Introducing null coordinates 
u = t — r and v = t + r for v > u, and compactifying them by U = arctanu and 
V = arctanu, we get 

iy= — i — ( -dV dU + \ sm\V - U)daA , on (-7r/2 < C/ < < 7r/2). 

cos^ V cos^ U \ 4 / 

Points at infinity with respect to the original coordinates have finite values with respect 
to the compactifying coordinates {V = tt/2 on U = — 7r/2), however, the physical metric 
fj in compactifying coordinates is singular at these points. This singular behavior can 
be compensated by a conformal rescaling with the conformal factor J7 = cos V cos U, 
so that the rescaled metric 

= n^r) = -dUdV + ^ sin^{V - U)da^ , 

is well defined on the domain {—tt/2 < U < V < tt/2) including points that are at 
infinity with respect to j]. We say that 77 can be extended beyond infinity. 

For Fig. 11.11 time and space coordinates t = {V + U)/2 and r = {V — U)/2 have 
been introduced. The resulting metric 77 = —dP + dr^ + sin^ r da^ is the standard 
metric on the Einstein cosmos R x S^. The embedding </> : R"* ^ R x 5*^ of the 
Minkowski spacetime into the Einstein cosmos is given by 

\ (^) + (^)) ' (^) " (^)) ' 

from M =M.^ = {t e (-oo,oo),f G [0,oo)} into M = R x = {t ^ (-00,00), 
r G [0,7r]}, suppressing the angular coordinates. The part of Ai that corresponds to 
the Minkowski spacetime is given by{|i + r| < tt ,\t — r\ <7r}. 

We see that the completion of is a manifold with boundary {t = ±{tt — r), 
r G [0,7r]} where the boundary points correspond to points at infinity with respect to 
the physical metric. Asymptotic behavior of fields on Ad can be studied using local 
differential geometry on this boundary where the conformal factor n = cos t + cos r 
vanishes. The part of the boundary without the points at r = 0, tt is denoted by 

= {t = ±(7r — r),r € (0, tt)}. The differential of the conformal factor does not 
vanish at J^, dD,\,y ^ 0, and consists of two parts J^^, each of them with the 
topology R X 5^. 

Fig. 11.11 shows curves that correspond to constant values of the coordinates t and 
f in the new representation. Each point of the diagram represents a sphere except 
the dashed vertical line segment connecting i"*" and which corresponds to the origin 



1.3 Conformal infinity 



13 




= const. 
= const. 



i 



Fig. 1.1: Penrose diagram of the Minkowski spacetime 

{f = 0}. Surfaces of f = const, approach timehke infinity = {t = ±7r,r = 0}, 
surfaces of f = const, approach spatial infinity i'^ — {t = Q,r ~ vr}. The figure is a 
Penrose diagram which is especially useful for depicting causal structures. In Penrose 
diagrams, radial light rays are straight line segments with 45 degrees to the horizontal. 
Null rays {U = const, and V — const, surfaces) reach y for an infinite value of the 
physical affine parameter along them, hence ^ is called null infinity. 

The idea is now to take certain properties of the asymptotic structure of the Min- 
kowski spacetime as being representative for the asymptotic structure of isolated sys- 
tems. We make the following assumptions about an asymptotically flat spacetime 
{A4,g): There exists a regular spacetime {A4,g), a sufficiently differentiable (say C^) 
function U on A4 and an embedding (f) : A4 ^ A4, which is conformal with conformal 
factor f2 such that g — {ilo (f))^'^(f)*g — cf)* {^l^"^ g) . The conformal factor satisfies f2 > 
on 4){AA) C Ai and the completion of 4){Ai) is a submanifold with non-empty bound- 
ary that is (90 (A^) = ^ 7^ 0, on which ri|.y — 0, 7^ 0. These assumptions 
are closely related to the notion of weak asymptotic simplicity |108j . In the following, 
we identify (/'(A^) with M. so that we write M. d M. or g — Q'^g. The triple {M,g, fl) 
is called conformal extension of {A4,g). Note that the conformal extension is not 
unique. 

The basic assumption here is that the rescaled metric g is regular across .y. A 
similar construction where this assumption is satisfied can be made for other explicit 
solutions that can be regarded as asymptotically fiat, for example the Schwarzschild- 
Kruskal solution (see Fig. 1 1.21 for the resulting Penrose diagram). However, we do not 
have a large class of explicit radiative solutions where we can test our requirements on 
the conformal extension. Therefore we do not know whether the above assumptions 
cover sufficiently general spacetimes so that we can take a conformal approach in the 
isolated system idealization of interesting astrophysical configurations . 

A strong support for the conformal compactification technique comes from the 
studies described in the previous section. Results achieved by coordinate dependent 
methods and limits to infinity can be derived very elegantly for spacetimes admitting 
a smooth conformal boundary using the conformal technique [107| I108j . 

To utilize the conformal compactification technique for calculating gravitational 
radiation we would like to answer: How general is the description proposed by Penrose? 
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The details of tlic conformal compactification technique have been motivated by 
expHcit examples or well studied assumptions about the fall-off behavior of certain 
fields as above and not by a detailed study of the full non- linear Einstein equations. The 
assumption on smoothness of J' is geometrical and includes only requirements on the 
conformal class \g\ . Smoothness of the rescaled metric g along the conformal boundary 

corresponds to a fall-off behavior for the physical metric g, which is however a 
solution to the Einstein equations that impose their own asymptotic behavior. It 
might be that the smoothness requirement on the conformal extension is too restrictive 
such that only a very special class of solutions satisfies it. Whether the geometrical 
requirement of smoothness of is compatible with the fall-off behavior of the metric 
fields as a solution to the Einstein equations is a delicate question. 

We would like to know by some general argument whether we have a large class of 
non-trivial, asymptotically flat, radiative spacetimes that admit a smooth conformal 
boundary such that we can apply the conformal compactification technique. To answer 
this question, the solution space to the Einstein equations needs to be studied with an 
emphasis on the asymptotic structure of gravitational fields. 

The available explicit solutions are not general enough to study the solution space 
by direct means, therefore we are led to abstract analysis. By sufhcient knowledge 
on properties of spacetimes, we can get some results in certain classes. We know, for 
example, that asymptotically flat, vacuum, stationary spacetimes admit an analytic 
compactification at null infinity [42,. For more general results that include radiative 
spacetimes, the initial value problem needs to be studied in a general setting. 

The analysis of the initial value problem for Einstein equations in the compactified 
picture is not straightforward, as the Einstein equations are not conformally invariant 
and compactification leads to formally singular equations at infinity. However, the 
equations are conformally regular as Friedrich showed by constructing a system which 
is equivalent to the Einstein equations for > and is regular for all values of the 
conformal factor so that the equations for the rescaled metric g can be analyzed on a 
conformal extension Ai including ^ where $7 = [13] . 

Friedrich analyzed the initial value problem for the conformally regular field equa- 
tions based on spacelike surfaces that extend smoothly through null infinity. These 
surfaces are called hyperboloidal as their asymptotic behavior is similar to the standard 
hyperboloids in Minkowski spacetime [54] . Let S be such a hyperboloidal surface that 
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Fig. 1.3: Penrose diagram for a hyperboloidal initial value problem. 

cuts J'^ in a spacelike two dimensional surface E = 5 fl J'^ and extends smoothly 
through ^+ as in Fig. 11.31 (the discussion of J'~ follows by time reversal). The anal- 
ysis of the hyperboloidal initial value problem revealed that if regular data for the 
conformally regular field equations is given on 5, the smoothness of J^"*" will be pre- 
served into the future of 5, at least for a while [54]. While this result is local in time, 
Friedrich also proved a semi-global stability result stating that, for small data, 
will admit future complete null generators and a regular timelike infinity |56j . The 
existence of regular data for the conformally regular field equations has been shown 
in [3]. The analysis of [3] showed further that, in general, even with the strongest 
smoothness requirements on free data, logarithmic singularities can arise at E unless 
some mild regularity conditions on the geometry of E are satisfied. 

These analytic results are important for practical purposes in numerical calcula- 
tions. One should be aware that the requirement of the existence of a smooth 
implies a restriction on the data. One may expect that this restriction does not ex- 
clude physically relevant solutions, but this is yet an unresolved question. Still, for a 
non-smooth some of the relevant structure may be recovered j351 1151j . It is an 
interesting open question whether one can deal with some mildly singular behavior 
numerically. We will study the hyperboloidal initial value problem in chapter [2l 

The analytic results just described show that the decision on smoothness of null 
infinity is made at spatial infinity Remarkably, our knowledge on the asymptotic 
behavior of solutions like Schwarzschild or Kerr can be applied to this question. By 
gluing techniques developed by Corvino and Schoen 36] , one knows that there exists a 
large class of non-trivial initial data which is Schwarzschild or Kerr in a neighborhood 
of spatial infinity. As stationary spacetimes admit an analytic compactification at null 
infinity, we know that the development of such data leads to regular hyperboloidal ini- 
tial data. Combined with Friedrich's results on the hyperboloidal initial value problem 
this leads to the existence of a large class of non-trivial, radiative, vacuum spacetimes 
with smooth null infinity [34]. These spacetimes, however, have a special asymptotic 
structure. The question still remains whether more general spacetimes exist. To an- 
swer it, the behavior of gravitational fields in a neighborhood of spatial infinity needs 
to be studied in detail. 

A detailed study of solutions in a neighborhood of spatial infinity is complicated 
by the fact that when the ADM-mass of the spacetime does not vanish point compact- 
ification at spatial infinity leads to a certain singular behavior. Physical fields do not 
admit smooth limits at the point . The limits depend on the spacelike direction along 
which one approaches j° [75l|4]. Therefore, to study fields at spatial infinity, it might 
be better not to represent it as a point. In [3 [TH I127j . a representation of spatial 
infinity as a unit timelike hyperboloid has been used. This representation of spatial 
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infinity, however, does not allow a study of the full system of the Einstein equations 
in a unified picture with spatial and null infinity. 

A regular finite initial value problem near spatial infinity has been formulated by 
Friedrich using the reduced general conformal field equations that he introduced in 
|58| . These equations are written in a conformal Gauss gauge based on conformal 
geodesies which are auto-parallel curves with respect to a Weyl connection [63l [7^ . 
In the conformal Gauss gauge, spatial infinity can be represented as a cylinder. The 
cylinder at infinity allows a detailed discussion of gravitational fields at spatial and 
null infinity in a single, regular setting that only depends on the conformal structure 
of the spacetime |61] |64j . 

Friedrich's analysis showed that in general, logarithmic singularities arise in a small 
neighborhood of spatial infinity. He obtained necessary regularity conditions on the 
Cotton tensor for these singularities to vanish. Later on, Valiente Kroon showed that 
these conditions are not sufficient and obtained further obstructions to smoothness 
of null infinity "141J . The question on necessary and sufficient conditions on Cauchy 
data for smoothness of is still open and a field of active research. A rich structure 
can be expected in a neighborhood of spatial infinity, as suggested for example by the 
existence of initial data whose development may lead to different smoothness properties 
on past and future null infinity [142| . 

For having numerical access to spatial and null infinity in a single setting, the 
reduced general conformal field equations are the only tool we have at present. A 
numerical implementation of this system needs to deal with certain difficulties. The 
regular finite initial value problem at spatial infinity distinguishes directions at spatial 
infinity. The regular representation of conformal data on a Cauchy surface crucially 
relies on a conformal rescaling of a distinguished tetrad with respect to which the fields 
are calculated. For numerical studies using this system, the geometry imposed by the 
cylinder needs to be implemented. We will discuss an implementation and tests in 
chapter [3l 

In this section, I described conformal techniques and some analytic results on their 
applicability to present the analytic foundation on which the numerical studies of this 
thesis are build. The studies are complementary to existing numerical methods but 
they are not an aim by themselves. A motivation to use conformal techniques in nu- 
merical calculations is given by the hope that interaction between mathematical and 
numerical methods will lead to fruitful studies of the solution space to the Einstein 
equations. Especially interesting seems the question, whether a mildly singular behav- 
ior of solutions at infinity can be treated by numerical methods so that conjectures 
on properties of solutions can be made which in turn might be studied by rigorous 
analyses. 

Another basic motivation is achieving efficient numerical simulations of isolated 
systems using controlled approximations so that reliable predictions on gravitational 
radiation can be made. Current numerical methods to calculate gravitational radiation 
have certain disadvantages that we will discuss in the next section to give an overview 
on the subject and to put the conformal approach in numerical relativity into its 
broader context. 
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1.4 The outer boundary in numerical relativity 

To study the development of large classes of data and to make quantitative predictions 
on highly dynamical, strong gravitational fields, we need to use numerical methods. 
The Einstein equations written in arbitrary coordinates are of no known type as they 
only determine the isometry class of the metric. To numerically calculate a solution 
metric with given initial data on some hypersurface, one typically reduces the Einstein 
equations to a hyperbolic system of partial differential equations. The reduced system 
is discretized and solved iteratively by some numerical algorithm. A reduction can be 
done in various ways which suggests a variety of methods in numerical relativity. 

Special properties of general relativity among other theories of physics, particu- 
larly the lack of a background and related difficulties like the non-linearity of the 
equations or the non-local feature of gravitational field energy, require a strong inter- 
play between mathematical and numerical methods not only on numerical analysis or 
analysis of partial differential equations, but also on differential geometry. The im- 
portance of abstract methods that can suggest well-defined techniques for numerical 
calculations has been recognized only recently. Today, the numerical stability problem 
for some astrophysically interesting classes of data seems to be solved. The numerical 
development of various binary black hole initial data, which was a big challenge for 
many years, can now be followed from the late inspiral through the coalescence phase 

There remain, however, many open problems that need to be studied. In this 
thesis we are mainly concerned with numerical studies of the asymptotic behavior of 
isolated systems. We will therefore concentrate on problems concerning the numer- 
ical treatment of the asymptotic region. To illustrate the computational domain in 
various methods that we discuss, the regions I and II in the Penrose diagram of the 
Schwarzschild-Kruskal spacetime from Fig. 11.21 will be plotted. We will also discuss 
the causal structure on the grid depicted by ingoing and outgoing null surfaces. We 
will assume that the interior of the event horizon is dealt with by the excision method. 
This corresponds to the introduction of a spacelike boundary inside the black hole so 
that no boundary data needs to be given. The requirement for the numerical treatment 
of such a spacelike inner boundary is numerical stability. There are other possibilities 
for treating the singularity that we ignore. 

1.4.1 Timelike artificial boundary 

The asymptotic region of isolated systems extends to infinity. To simulate such space- 
times numerically on a finite grid, the most common method truncates the computa- 
tional domain by introducing an artificial, timelike outer boundary which introduces 
certain problems that we discuss in the following. Fig. 11.41 shows a typical foliation 
of the Schwarzschild-Kruskal spacetime along with timelike surfaces that may act as 
artificial boundaries. The diagram has been calculated using the ingoing Eddington- 
Finkelstein coordinates described in Appendix [Al 

The introduction of an artificial, timelike, outer boundary implies an initial bound- 
ary value problem for a hyperbolic reduction of the Einstein equations which should 
be well-posed. The requirement of well-posedness is not just a mathematical subtlety 
but is also physically motivated. It means that unique solutions are required to exist 
that depend continuously on initial data. This is strongly related to the requirement 



1.4 The outer boundary in numerical relativity 



18 



outgoing 

i ingoing 

Ts = const. 

rs =0 i+ I;-,; I ^ . .J V X 7^ 
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of predictability on a physical theory. 

Further, the solutions to the reduced system should correspond to solutions to the 
full system in the sense that the vanishing of the constraints on an initial hypersurface 
is preserved into the future of that surface. Therefore, the boundary conditions at the 
outer boundary should preserve the vanishing of the constraints. Such a well-posed, 
constraint preserving initial boundary value problem for Einstein equations has been 
constructed by Friedrich and Nagy [68] ■ Their frame-based, first order formulation 
allows the adaptation of the evolution system to a timelike outer boundary such that 
the constraints are propagated by maximally dissipative boundary conditions without 
the need to restrict the boundary data. However, their reduction has not yet been used 
in numerical simulations. Currently, there are many studies on the the initial boundary 
value problem for the Einstein equations using methods that are simpler to implement 
in numerical calculations but do not solve the full problem P [71 [88 l [9T t fTTT l fTTSl , [121] . 

Another difficulty concerns the choice of boundary data. We see in Fig. 11.51 that 
ingoing characteristics enter the spacetime from the outer boundary. One might wish 
to prescribe data on these characteristics that would correspond to a no incoming 
radiation condition. However, the timelike outer boundary surface is geometrically 
arbitrary, without invariant meaning. Apart from the problem related to its gauge 
dependent choice, it seems difficult to give a precise definition of gravitational radiation 
on such a surface. This difficulty is related to the lack of a quasi-local energy concept in 
general relativity. Besides, such a no incoming radiation condition will not in general 
correspond to the transparency of the artificial boundary due to backscattering of 
gravitational radiation, which is related to the non-linearity of the theory. To give an 
example for potential problems, we mention that it has been proven that certain types 
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of commonly used boundary conditions in numerical relativity destroy the tail behavior 
[10]. Nonetheless, there are detailed analyses on how to deal with these problems in 
the linearized regime [551 127] • 

The problems described above can be treated numerically provided that the outer 
boundary is sufficiently far away so that linearized analysis is a good approximation. 
When the outer boundary is far enough, a numerically stable treatment might be 
good enough as one can trust the solution up to numerical accuracy on the domain 
of dependence of the initial hypersurface. As can be seen in Fig. 11.41 the domain of 
dependence grows when the outer boundary is put farther away. A difficulty arises 
due to the slow 1 jr fall-off rate of the radiation field. To gain a factor of 2 in accuracy, 
the outer boundary needs to be pushed to twice the distance which requires 2^ — 
8 times the number of grid points on a homogeneous grid in a simulation in three 
space dimensions. Therefore, today's numerical codes use inhomogeneous grids. By 
using mesh refinement techniques one can put the outer boundary to as far as 700m 
[9l[97] . Mesh boundaries, however, cause numerical errors and there is still a problem of 
efficiency on foliations that approach spatial infinity with increasing physical distance. 
Such foliations waste computational resources, as one needs long time evolutions to 
calculate the emission of radiation in the far-field zone. 

This brings us to the radiation extraction problem. In the analysis of numerically 
generated spacetimes, one calculates certain quantities interpreted as gravitational ra- 
diation along a timelike surface representing a family of observers far away from the 
source whereas is the natural place to measure emission [50, . As there is no unam- 
biguous definition of gravitational radiation at a finite distance away from the source 
one uses approximative methods. A common method, called the Regge-Wheeler-Zerilli 
method assumes that the full metric in the extraction region is a perturbation of the 
Schwarzschild spacetime in a certain gauge and tries to extract radiation information 
in a coordinate invariant way |115l I156j . Another method relies on calculating Weyl 
scalars with respect to a special tetrad class, called the quasi-Kinnersley frame, that 
can be constructed in an invariant way in spacetimes of Petrov type D [99j . These 
methods rely on linearized approximations. As shown in [106) . different extraction 
methods can deliver different results in an accurate code so that numerical error is 
dominated by systematic error from the extraction method. Note that this error can 
not be estimated by convergence tests. Even in the continuum limit of infinite resolu- 
tion, the calculation of gravitational radiation will have the same systematic error. 

To estimate the extraction error, one can put the extraction surface farther away 
so that the systematic error decreases typically with 1/f. However, there are limita- 
tions on the choice of the observer location that also limit the accuracy of radiation 
extraction. On the one hand, one should not put the observer too close to the sources 
because the assumptions underlying the extraction method are not valid in the strong 
field regime. On the other hand, one should not put the observer too far, because 
then the contamination from the outer boundary hits the extraction surface before 
the waves reach the observer. A further difficulty arising from putting the extraction 
surface far away is due to numerical dissipation which lowers the accuracy of the ex- 
traction as the waves propagate slowly along the grid. Today, the observer location 
for wave extraction is typically set to about 50m [l6l [97] . 

Although recent years saw quite impressive advances in the problems mentioned, 
one can say that there is still a lack of appreciation of global or semi-global considera- 
tions. While the problems we alluded to can be regarded as subtleties and not crucial 
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for the detection of gravitational waves, the approximation errors might turn out to be 
relevant for gravitational wave astronomy where accurate waveforms are needed that 
might in principle allow us to decide on different equations of state for inspiralling 
neutron stars or to extract the remnant gravitational radiation from the big bang. 
Such questions require accurate predictions under controllable approximations. As 
large parameter studies need to be done, the efficiency of the codes will also play an 
important role. 

Beside the requirement of an unambiguous radiation extraction we can deduce the 
following requirements on a smooth foliation for a clean boundary treatment just by 
considering the causal structure on the grid as depicted in Fig. 11.51 

• No characteristics should enter the computational domain from the boundary. 

• Outgoing characteristics should leave the computational domain through the 
boundary. 

Note that both requirements are fulfilled for excision along a spacelike inner boundary 
inside the black hole. 

In the following we will review alternative methods that have been suggested to 
deal with the problems of the artificial timelike outer boundary approach. 

1.4.2 Coordinate compactification at spatial infinity 

To avoid problems related to the treatment of the outer boundary within the standard 
3+1 approach, one can use a compactifying coordinate system on the Cauchy hyper- 
surfaces. Then the outer boundary of the computational domain is at spatial infinity 
where natural boundary conditions can be used [33l 1113] . 

To elucidate the idea, take an asymptotically flat, two dimensional spacetime metric 
and a time function t such that t = const, surfaces are Cauchy surfaces. Introduce 
on these surfaces a coordinate r. We write the metric on the coordinate domain 
{t G (— oo, oo), f S [0, oo)} as 

g — cjtt dt^ + 2gtf- dt df + gff df^ . 

To avoid a timelike outer boundary, we can introduce a compactifying coordinate r by 
setting r{r) = r/(l + f). The coordinate transformed metric takes the form 

g ^ gtt dt'^ + 2gtr dt dr + cjrr dr'^ , where = — g^^^—Ji:i—, (1.6) 

(1 — r)^ (1 — rj* 

The coordinate value {r — 1} corresponds to spatial infinity i° as the surfaces t = 
const, are Cauchy surfaces. While the transformed metric components gtr and grr are 
singular at this point, the untransformed ones attain their Minkowski values for an 
asymptotically flat spacetime. A numerical code can evolve the untransformed com- 
ponents by substituting the relation p.6|) into the Einstein equations. The equations 
become formally singular at {r = 1} where the following outer boundary conditions 
can be applied 

~gtt{t,l) = -I, gtf{t,l)^0. ~gff{t,l) = l. 

The outlined method solves the Einstein equations in a domain that is not bounded 
by a timelike dashed curve in region I of Fig. [Ll but by the point i°. Note that 
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Fig. 1.6: Causal structure on the 
grid bounded by spatial infinity and a 
spacelike surface in the interior of the 
event horizon. 



Fig. 1.7: Causal structure on the same 
grid at a later time. 



although is a point, in the coordinate compactification at spatial infinity it is blown 
up in terms of the coordinates as illustrated in Fig. 11.61 and Fig. 11.71 In [1131 1114j 
spatial infinity becomes a cylinder with cubical cuts due to the use of a Cartesian 
coordinate system. 

While this method avoids the introduction of an artificial boundary, it does not 
solve the problem of radiation extraction. Using compactifying coordinates on Cauchy 
surfaces is very awkward for studying radiation. The grid velocity of outgoing char- 
acteristics decreases to zero towards the outer grid boundary. In simulations, this is 
seen by a slowing down of outgoing waves. Gravitational radiation, traveling along 
outgoing null rays in Fig. 11.61 and Fig. 11.71 can not leave the numerical domain and 
piles up during the evolution near the boundary. This leads to instabilities as the 
waves can not be resolved after some time. Dissipation can be used to deal with such 
instabilities, but this leads to a loss of accuracy. Radiation extraction is done at some 
timelike, gauge dependent surface at a finite distance away from the source as in the 
artificial timelike boundary case which implies similar problems. Extraction at outer 
domains is not accurate due to the decrease of wave resolution at large spatial dis- 
tances (see |113[ 1118] for detailed discussions). In addition, errors generated in the 
vicinity of the outer boundary can travel along ingoing null rays depicted in Fig. 11.71 
and contaminate the interior solution. 

In this method, our first requirement that no incoming characteristics should enter 
the computational domain is fulfilled, while the second requirement is not. The use 
of hyperboloidal foliations instead of Cauchy foliations might lead to a method that 
satisfies the second requirement too. To my knowledge this idea has not yet been 
tried out. It is an interesting open question whether the coordinate compactification 
technique can be applied using hyperboloidal foliations. 
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Fig. 1.8: Penrose diagram of the CCM- Fig. 1.9: Causal structure on the grid 
approach with a matching boundary at for CCM with a matching boundary at 
fs = 3m. fs = ^ITT- 



1.4.3 Characteristic evolution and matching 

The characteristic approach is based on studies by Bondi, Sachs and others [2Tlll20| . 
For the treatment of the outer boundary in numerical calculations with the charac- 
teristic approach, one uses compactifying coordinates along outgoing characteristic 
surfaces with respect to which ^+ is at a known coordinate location and rescales the 
metric with an appropriate conformal factor 134]. In a certain conformal gauge, this 
procedure leads to equations which are formally singular but have a regular limit at 
that can be calculated by numerical methods. 

The number of variables in a characteristic evolution scheme reduce dramatically 
compared to the corresponding version of a Cauchy problem because the equations 
do not involve second time derivatives. Also the structure of the equations is simpler 
and there are no elliptic constraints so that initial data are free [152]. These features 
make the characteristic approach very attractive. This approach has been successful in 
cases where null foliations smoothly cover the numerical domain, such as a single black 
hole or a relativistic star. Unfortunately the coordinates are not flexible and there is 
little gauge freedom one can use to avoid coordinate singularities. The main problem 
with this approach is the fact that characteristic foliations are not well-behaved in 
regions of strong dynamical gravitational fields due to formation of caustics in the 
light rays generating the null hypersurfaces [71]. Therefore, a modification has been 
suggested which matches a calculation in the interior based on a spacelike foliation 
to a characteristic calculation in the asymptotic region where the fields become weak 
and a null foliation can be expected to smoothly cover the simulation domain. This 
approach is called Cauchy-characteristic matching (CCM) [191 1152] . 

The conformal diagram for a CCM evolution scheme has been plotted in Fig. 11.81 
Note that the solid curves represent pieces of spacelike Cauchy surfaces as in Fig. 11.41 
Outgoing characteristics from the interior calculation are developed up to ^+ and 
ingoing characteristics to the interior are calculated by interpolating the solution in 
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the characteristic region to a timelike curve. This timelike curve represents an outer 
boundary for the Cauchy evolution and an inner boundary for the characteristic evo- 
lution. The causal structure on the grid resulting from the matching has been plotted 
in Fig. \TM 

Both the inner and the outer boundaries fulfill our requirements for a good nu- 
merical boundary treatment, however, the foliation is not smooth. The main difficulty 
with CCM is the matching along the timelike boundary between the interior and the 
asymptotic region, where the causal nature of the foliation changes. Stability prob- 
lems caused by interpolation between a Cartesian code solving a Cauchy problem and a 
spherical code solving a characteristic problem using different sets of variables impeded 
a further development of this approach. One can hope that with recent improvements 
on the treatment of timelike boundaries, it will be possible to overcome the difficulties 
the CCM approach has been facing in the past. 

In the next subsection we will see that the required behavior at both boundaries 
can be achieved smoothly without changing the causal nature of the foliation or the 
variables of the reduction. 

1.4.4 Conformally regular field equations 

The conformally regular approach in numerical relativity is based on analytic studies 
by Friedrich |531 1541 156j and it started with numerical studies in spherical symmetry 
by Hiibner |79| . In this approach one solves numerically a hyperboloidal initial value 
problem for the conformally regular field equations on a domain illustrated in Fig. 11.101 
(see [51] [82] for reviews) . The use of compactified hyperboloidal foliations is promising 
because they combine advantages of Cauchy and characteristic foliations. On the one 
hand, instead of approaching spatial infinity as Cauchy surfaces do, these surfaces 
reach null infinity which makes them suitable for unambiguous radiation extraction, 
and on the other hand, in contrast to characteristic foliations, hyperboloidal foliations 
are spacelike everywhere and therefore as flexible as Cauchy surfaces. 

An exemplary causal structure on the grid is depicted in Fig. 11.111 for constant 
mean curvature slices in the Schwarzschild-Kruskal spacetime in "^-fixing coordi- 
nates, i.e. in coordinates in which the spatial coordinate location of has been fixed 
as described in 12.2.21 Both requirements for a good boundary treatment arc fulfilled. 
Outgoing characteristics leave the computational domain through null infinity where 
a rigorous definition of gravitational radiation allows us to construct an unambigu- 
ous numerical radiation extraction method and there are no incoming characteristics 
so that no boundary conditions are needed. The foliation is smooth throughout the 
simulation domain and very flexible. 

A difficulty with the conformally regular approach in numerical calculations is 
that the equations include, among others, evolution equations for the conformal fac- 
tor which results in a solution dependent representation of ^+ = {$1 = 0}+. Note 
that in Fig. II. 11) ^+ has been fixed to a spatial coordinate location. In the metric 
based conformally regular field equations, however, the location of J^"*" on the grid is 
not known a priori. The numerical boundary does not coincide with the conformal 
boundary ^+ which leads to problems of efficiency and requires a numerical bound- 
ary treatment outside the physical spacetime M [81]. A gauge condition to fix 
to a spatial coordinate location has been suggested by Frauendiener in [49] in the 
context of frame-based conformally regular field equations. To my knowledge, this 
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Fig. 1.10: Penrose diagram of a hyper- 
boloidal initial value problem on the extended 
Schwarzschild spacetime. 



Fig. 1.11: The corresponding 
causal structure on the grid. 



gauge has not been implemented in a general scheme in which corresponds to the 
outer grid boundary due to technical complications with numerical implementation of 
a frame-based evolution system that requires spherical grid topology. 

Another difficulty is that the conformally regular field equations are significantly 
larger than usual formulations of Einstein equations including many additional evolu- 
tion variables subject to constraints. Due to the large number of constraint equations, 
numerical errors require a stronger control on constraint propagation properties of the 
system. As there is not enough numerical experience with the equations, one can not 
use established methods easily to deal with the encountered instabilities. While these 
problems are not of a principal nature, they have made progress in the conformally 
regular approach difficult. 

Instead of trying to overcome the described difficulties within the context of the 
outlined methods, we will discuss in the following chapters two novel conformal ap- 
proaches for the numerical calculation of asymptotically flat spacetimes. 

The approach taken in chapter [3] falls into the category of the conformally regular 
approach. It is based on the system of reduced general conformal field equations which 
has certain advantages over the systems that have been tried in the conformally regular 
approach until now. Using this system, we will numerically study a neighborhood of 
spatial infinity including a smooth piece of null infinity. 

Once a full neighborhood of spatial infinity has been calculated, including a piece of 
null infinity, one can evolve the system further by solving a hyperboloidal initial value 
problem. For accurate numerical calculations of gravitational radiation along null 
infinity, however, the reduced general conformal field equations do not seem to be the 
appropriate tool. The reason is that the underlying gauge based on timelike conformal 
geodesies leads to a loss of resolution in the physical domain. A further motivation 
to construct a different method for the numerical solution of a hyperboloidal initial 
value problem is the wish to employ the extensive experience in numerical relativity 
gathered over the last decades. In chapter [21 I will present a method in which null 
infinity can be included in the computational domain for a common reduction of the 
Einstein equations. In this method the location of null infinity is fixed on the numerical 
grid so that no resolution loss appears. 



Chapter 2 

Null Infinity 



The discussion on numerical outer boundary treatments in the introduction suggests 
that one should include null infinity in the computational domain for a proper treat- 
ment of the outer boundary as well as for unambiguous radiation extraction. In this 
chapter, we will see how J"^ can be included in the computational domain for a 
common reduction of the Einstein equations. 

We start with a study of a certain class of spherically symmetric hyperboloidal 
surfaces in the Minkowski spacetime to develop an intuition about their behavior. 
We observe in spherical symmetry that a conformal factor together with a radial 
spatial coordinate r can be chosen such that is given by {r = 1} as plotted in 
Fig. 11.111 Then we address the full problem without symmetry assumptions. We 
see for a common reduction of the Einstein equations how .^"'"-fixing can be achieved 
by prescribing the representation of a preferred conformal factor in terms of suitably 
specified coordinates in a hyperboloidal initial value problem. We present numerical 
tests of the method in spherical symmetry. The chapter ends with a discussion. 



2.1 Hyperboloidal surfaces in Minkowski spacetime 

A hypersurface iS in a spacetime (^M. , g) can be given as a level set of a scalar func- 
tion <I>(p) = const, which satisfies ci$ 7^ for all p £ S. We assume {M,g) to be 
time-orientable, with metric signature (— , +, +, +) and coordinates {x^^}. We set the 
positive a;°-direction to be the future direction. The unit normal to the hypersurface 
is given by 

""''■"WW.- 

We require the hypersurface to be space-like so that the unit normal is time-like, 
g^iyU^n}' = — 1, and we choose to be future-pointing, > 0. If the level sets of $(a;) 
define a local foliation in A4, we can introduce a new time coordinate t{x) — $(a;) and 
construct a new coordinate system on a certain domain by choosing space coordinates 
on one of the hypersurfaces 5$ = {$ = const.}. We define the lapse function a and 
the shift vector f3^ in this coordinate system by the decomposition 9f = an^ + (3^, 
with g{n, (3) ~ 0. The induced Riemannian metric on St is given by h^^, = g^^, + 11^^11^. 
The extrinsic curvature K^^, given by K^^ = h^W\n^, is a measure for the variation 
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of the unit normal along the hypersurfaces St ■ The mean extrinsic curvature K is the 
trace of the extrinsic curvature, 



1 



-9 

In our notation positive K means expansion. The mean curvature of a hyperboloidal 
surface S approaches a strictly positive value at . We will discuss some examples of 
spherically symmetric hyperboloidal surfaces in the Minkowski spacetime (M^jiy) with 
the metric fj given as in (jl.5p . The mean curvature K of such a hypersurface with unit 
normal h'^it^f) becomes 

K = din* + dfff + ln^. (2.2) 



2.1.1 Standard hyperboloids 

The standard hyperboloids in the Minkowski spacetime can be given as level sets of 
the function 

$(t,f)=f-r2. (2.3) 

We investigate the surfaces 5$ = {<& = const.} strictly inside the upper light cone, 
i > f, where they are spacelike and $ > 0. The unit future-directed normal (|2.ip and 
the mean curvature (j2.2p on 5$ read 



Standard hyperboloids given by different values of $ have spatially constant mean 
curvature (CMC), however, if we introduce $ as a foliation parameter we do not get 
a CMC foliation as the value of K varies from slice to slice. 

The embedding of the Minkowski spacetime into the Einstein universe from 11.31 
allows us to study the asymptotic behavior of standard hyperboloids by using local 
differential geometry. We get 

$ {i{V, U),r{V, U)) = tan V tan U. 

The intersection I]$ = 5$ n is independent of <I>, 

U\^^ = arctan(<i>coty)L,_„ = 0, 

as seen in Fig 12.11 Therefore $ is unsuitable as a foliation parameter for J^~^ . It 
determines the angle a of the cut at S = {U = 0,V = 7r/2}, 



tanJJ tanV^ ^ dU 

cos^ V cos-' U dV 



tan U cos^ U cos^ U 
= = — <j) 

tan V cos^ V sin^ V ' 



dU 

tana = 

dV 



$ a = arctan$ = arctan ( ) . (2.4) 



We can use this relationship to find the hyperboloid that intersects ^+ in a certain 
angle. The hyperboloid which is "parallel" to the t — surface in the above embedding 



2.1 Hyperboloidal surfaces in Minkowski spacetime 



27 




Fig. 2.1: Standard hypcrboloids. 




V 



Fig. 2.2: Cutting angle at J^+. 



has a cutting angle a = 7r/4 which implies $ = 1. The surface 5i = {$ = 1} has been 
plotted in Fig. [O and Fig.O 

We see from (12.41) that for large values of if, the angle approaches 0. This suggests 
that K can be used as an intuitive measure for how close the surfaces are to null 
surfaces. In subsequent sections we will see how this aspect plays a role in numerical 
calculations. 



2.1.2 Simple hyperboloidal foliations 

A hyperboloidal foliation can be constructed by translating a standard hyperboloid, 
given by = t — f ^ with a S M, along the timelike Killing vector field di of the 
Minkowski spacetime. The translated surfaces satisfy a? — [i— — f^. So, they are 
given by level sets of 

= Va^ + (2.5) 
The unit future-directed normal and the mean curvature read 



f 2 a- 



a 



In this case, each surface not only has a spatially constant mean curvature, but also 
different surfaces of the foliation share the same K. The value of the mean extrinsic 
curvature does not depend on $ but on a real parameter a that can be fixed freely. 
We embed the surfaces into the Einstein universe as before to get 



<S>{V,U) = i(tany + tanC/) - \j + ^(tan V" - tan C/)2. 

To study their asymptotic behavior near ^+ where V 11/2 and cot ^ we write 

2 tan [/ - 2 cot 
U) = , =■ 

1 + tan [/ cot y + a/ 4a2 cot^ V + (1 - tan C/ cot F )^ 
We get a Taylor expansion in cot V near ^+ 

^{V,U) = i&nU - coiV - a^i&nU co\?V + 0{coi^V) for F ^ 7r/2. 
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Fig. 2.3: CMC foliations of the Minkowski spaeetime with K = {15, 3, 1.5}. 



The cut at J^"*" depends on $ only 



$1 



— tan U 



U 



— arctan$. 



The angle of the cut can be calculated by 



tana 



dU 
dV 



a = arctan 



if2(l + $2) 



(2.6) 



Fig. l2.3l shows three CMC foliations in the Penrose diagram of the Minkowski spacetime 
with the same set of values of $, but different values of K. Each plotted surface is 
spacelike, so their angle is smaller then 45 degrees to the horizontal. For small K 
some of the surfaces look similar to Cauchy surfaces while for large K the surfaces 
rather look like null surfaces. The diagrams illustrate the sense in which we say 
that hyperboloidal surfaces mediate between Cauchy and characteristic surfaces. This 
behavior affects also the dynamics of wave-like solutions. The propagation speed of 
waves on the grid in adapted coordinates depends among others on the choice of K as 
can be seen in Fig. 12.41 It is a nice feature that by choosing the mean curvature we 
can control the behavior of the surfaces in the interior without changing their cut at 

For a given foliation with some K, we can find the surface which is "parallel" 
to the Cauchy surface t = with the slope a = 7r/4 at From ()2.6|) we get 



±\/a^ — 1 = iyQ/if^ — 1. Such surfaces exist only for foliations with K < 3 
which can also be seen in Fig. 12.31 

Another interesting set of foliations can be constructed by requiring that each 
surface of the foliation intersects ^+ with the same angle. Choosing a = 7r/4 we 
substitute = 1 + (f)^ into ((231) to get 



$(f, r) 



2i ( 



The mean curvature reads by construction K = 3/vr+$^. As in the case of the 
standard hyperboloids, each surface has a spatially constant mean curvature but K 
depends on the foliation parameter <i>. In this case we get a foliation. The embedding 
results in 

$ = - cot(y + [/) = - cot(2i), 

where we used the compactifying time coordinate t — {V + U)/2. As might have been 
expected, we get the t — const, surfaces. 
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2.2 c^^-fixing in spherical symmetry 

As seen in Penrose diagrams like Fig. l2.3i is an ingoing null surface. This property 
can be made manifest in suitable coordinates such that the null generators of 
converge. Such coordinates can be useful in numerical studies for discussing global 
properties of spacetimes, for example with respect to the existence of a regular point 
corresponding to timelike infinity i"*" |80| . These coordinates, however, do not seem 
to be convenient for calculating gravitational radiation accurately as they lead to a 
loss of resolution in the physical part of the conformal extension. Besides, when the 
outer boundary is a timelike surface in the unphysical spacetime, our requirements for 
a boundary treatment from 11.41 are not fulfilled. 

It has been suggested [31 HHl IHll 1155] that conformal compactifications in which 
is kept at a fixed spatial coordinate location might be suitable for numerical 
calculations. This is also suggested by illustrations of the causal structure as in Fig. II. 91 
or Fig. 11.111 In this section, we discuss the general construction of such ^+-fixing 
coordinates for spherically symmetric hyperboloidal foliations. The explicit examples 
concentrate on CMC foliations. 

We introduce some conformal transformation rules. Assume a time coordinate t 
and space coordinates {a;'''}^=i^2,3 have been chosen in M.. We write the conformal 
rescaling in terms of variables of a 3+1 decomposition in the given coordinate system 

g = {-o? + h^sfi^P^) dt^ + 2h^sP^ dt dx^ + h^s dx''dx^ = fl^g = 

= (^-{naf + n%sP''P^) dt^ + 2n%sP'' dtdx^ + n%sdx''dx\ 

and deduce the transformation rules 

The unit normal — (9f — P^) transforms as = jj"-''- ^^'^ transformation of 
the covariant derivative under the conformal rescaling g — fl^g reads (144j 

V + s{Q-'dn), s^jn-^dn) = ^{d^^n, + s^n^-g^pg^Mx), (2.7) 

where fl^ := (dfi)^ = V^fi. By the above transformation rules and by the definition 
of the extrinsic curvature -FT^^ = h^V\n^, we derive K^^i, — r2(ii'^,^ + ri„/7,^^), where 
r2„ = n°'Vla- The mean curvature transforms a,s K = Q.^^{K + 3ri„). We note that a 
CMC foliation with respect to the physical metric g does not necessarily imply a CMC 
foliation with respect to the rescaled metric g. 

The transformation rules above have been written in a given coordinate system. 
For ^+-fixing, we need to construct the coordinate system together with the choice 
of the conformal rescaling. This can be done explicitly in spherical symmetry. The 
physical line element in spherical symmetry in coordinates {i, f, x^} can be written as 

g = (-52 + /i2(^f)2^ dt^ + 2hji3^dtdf + hjdf^+f^da\ (2.8) 

The lapse a, the shift component /?'', and the spatial metric function hf are func- 
tions of the coordinates {t, r) only. We assume that the metric p.8p admits a smooth 
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conformal boundary and that the time coordinate t is such that t = const, hyper- 
surfaces are hyperboloidal. We do not compactify the time direction. The conformal 
compactification, g — fl^g, can be done such that 

n'^ihjdr^ + f^da^) = hi dr^ + dcr^ (2.9) 

with respect to a compactifying radial coordinate r. Note that we have some freedom 
here. One can require for example hr — 1 which leads to r being the proper distance, 
but then the conformal factor is determined by the above relation and the radial 
coordinate transformation can not, in general, be written in explicit form. By keeping 
hr, we have the freedom to prescribe the representation of a conformal factor in terms 
of a suitable compactifying radial coordinate r and the coordinate transformation 
can be made explicit. The relation (|2.9p implies for a given representation il(r) a 
coordinate transformation f — r so that dr — {il — r dr. Then the spatial 

metric function transforms a.s hr = hr = (fl — r n')fl~^hf. For regularity of this 
conformal compactification, hf{t, f) needs to have a specific asymptotic fall-off behavior 
for f — > CX3 on the hyperboloidal surfaces of constant t. 

A convenient representation for the conformal factor is f2 = (1 — r). This is not a 
good choice at the origin, but we are interested in the asymptotic region. A possible 
representation of a conformal factor in the interior will be discussed in 12.4.41 in con- 
nection with a hyperboloidal initial value problem. Our choice fixes the compactifying 
coordinate r via 



r 



dr dr 



which implies 

a = Via, = n^fi^ hr = n-^hf. (2.11) 

In the following we will explicitly construct a c^+-fixing conformal compactifica- 
tion for Minkowski and extended Schwarzschild spacetimes using a CMC foliation. 
There are many other possibilities for constructing simple hyperboloidal foliations, for 
example constant scalar curvature surfaces |104j . 

2.2.1 The Minkowski spacetime 

We use the CMC foliation (|2.5|) to introduce a new time coordinate t = ^{t, r) for the 
Minkowski spacetime (M^, fj). We have di = dt + r/ \/a^ + W dr, so that the Minkowski 
metric ()1.5p becomes 

77 -dt"^ "^^ dtdf + — - dr^ + da"^. 

The variables of the 3-1-1 decomposition read 

a = —y/a^ + r'^, — — — a, hf — —. 
a a a 

Conformal compactification, r/ = VP'fj, with = 1 — r using (|2.1ip results in 

1^ = -n^dt' - dtdr + f ^ dr^ + r^da^. (2.12) 
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Fig. 2.4: Causal structures on the grid for CMC foliations of Minkowski spacetime 
with K = {15,3, 1.5}. 



The variables of the 3+1 decomposition are 
a{r) = — V fl"^ 



r2 = n< 



^ a = ^^'/3^ 

a 



a 



With our choice of the conformal factor we get ri„ — (3"^ /a — —r/a = —Kr/3. Using 
the transformation rule of the mean curvature under conformal rescalings we derive 
K = il^^{K + Srin) = fl^^{K — Kr) = K. In this case the mean extrinsic curvatures 
with respect to the physical and the rescaled metric are equal so that a CMC foliation 
of the physical Minkowski spacetime gives a CMC foliation of the conformal extension 
that we constructed. 

The causal structure for rj plotted in Fig. 12.41 is equivalent to the causal structure 
for fj. The radial compactification allows us to include ^+ in the grid domain and the 
rescaling of the metric results in a regular conformal extension. We see that the value 
of K has a strong effect on the grid speed of outgoing characteristics. 



2.2.2 The extended Schwarzschild spacetime 

To find a representation of the Schwarzschild spacetime based on hyperboloidal sur- 
faces, we use the family of spherically symmetric CMC surfaces in the extended 
Schwarzschild spacetime constructed in [96] (see Appendix IA.1.21 for the definition 
of the surfaces and [73] for a numerical study). In coordinates adapted to the CMC 
slicing, the standard Schwarzschild metric is obtained in the form 

/ 2m\ ,3 2(J(f)-C) , o , , 

0, = - 1 dt^ - --XJ- !- dtdr + ^ df^ + da^, 

\ r J p{f) p2(^) 

where C is a constant parameter and we have defined 

J(f) := I r\ Pir) ^( J(f) - C)^ + (^1 - 

The constant mean curvature of the surface St = {t = const.} is denoted by K. The 
unit normal to the hypersurfaces St is 

P{r) r 
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Fig. 2.5: Causal structures on the grid for CMC foliations of the Schwarzschild space- 
time with C = 2, m = 1 and K = {1,0.3,0.07}. See Fig. 12.71 for the corresponding 
Penrose diagrams. 



The variables of the 3+1 decomposition read 

. P{f) ~, J{r) - C , ~ 1 

a ~ — 7y—, p = ^ a, hf. — —. 

f f a 

Conformal compactification, gs = ri^^s, as in (|2.9p with Vl = 1 — r and (|2.10p 
results in 

= - 1 n^dt^ \ , dtdr + dr^ + r^da^. (2.13 

\ r J Pyf) P^[r) 

where J(r) = J{f{r)) and P{r) il^ P(f(r)), or 

J(r) = I r^ P{r) - Cn'f + (^l - ^^^^^^) (1 - r)\'^ ' . 

The variables of the 3+1 decomposition read 

a=—^, = a, hr = -. (2.14) 



The mean curvature reads 

K = K ^ TT- 



The slicing of the rescaled metric (I2.13P is not a CMC slicing for C ^ 0, but the values 
of K and K are equal at 

Different choices of K by constant C have a similar effect on the grid speed of 
outgoing characteristics as in the case of CMC slicings of the Minkowski spacetime. 
Large values of K correspond to high grid velocities of outgoing characteristics as 
seen in Fig. 12.51 For low values of K, the causal structure on the grid resembles 
the causal structure on Cauchy surfaces (Fig. II. 7p in that the grid speed of outgoing 
characteristics decreases. In contrast to a Cauchy foliation, this speed stays larger 
than zero allowing outgoing characteristics to leave the spacetime through the outer 
grid boundary. 
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2.3 Hyperboloidal evolution with a prescribed Q > 

We have seen in spherical symmetry that we can prescribe the representation of a 
conformal factor in terms of a suitably chosen radial coordinate which can be deter- 
mined to achieve explicit ^+-fixing. The freedom to prescribe the representation of a 
conformal factor might be a special feature of spherical symmetry. We will see in the 
following, however, that one can construct general solutions to the Einstein equations 
where the representation of a conformal factor can be prescribed freely in terms of 
suitably specified coordinates. 

Consider the transformation behavior of the Ricci tensor under conformal rescalings 
g = VL^g with f7 > [m] 

= - ^ (2 V^V^l^ + nVLg^^) + ^(Va^^) V^f75^^. 

The Einstein tensor, G^^ = i?^,^ — ^g^^R, transforms as 

G^M = Gf,M - ^ (V^V.r! - U^g^,,) - ^{V g,,,. (2.15) 

As seen by this relation, the Einstein vacuum field equations for the physical metric, 
Giw[g] = 0, are equivalent to a similar system 

G^,M = T^.[r!] := (V^V^f! - Ung^,) - A(V;,r!)V^l] g^., (2.16) 

for the conformally rescaled metric g. We want to formulate a well-posed initial value 
problem for this system. 

As mentioned in [62] there are two difficulties when we try to solve (|2.16p for 
a conformally compactified metric g. The first one is caused by the terms involving 
divisions by factors of f2 which are formally singular at . These terms attain regular 
limits at J^^ if the spacetime admits a smooth conformal boundary. In the following 
we will assume that hyperboloidal initial data has been chosen whose development 
admits a smooth conformal boundary as indicated in 11.31 The question how these 
limits can be calculated numerically will be discussed in l2.5l 

The second difficulty is related to the question of how to determine Vt. The con- 
formal extension is what we are solving the equations for and the conformal factor is 
related to the asymptotic structure of the solution metric, therefore it must be deter- 
mined jointly with the metric. Although it seems difficult to control the evolution of 
the conformal factor we will see that the conformal invariance of our evolution system 
(|2.16p combined with certain geometric properties of allow us to prescribe the 
conformal and the coordinate gauge in a way that is kept at a spatial grid coor- 
dinate, the representation of a suitable conformal factor in terms of grid coordinates 
is known a priori and each of the formally singular terms in (|2.16|) attains a regular 
limit at . 

In this section we concentrate on the case of a positive, sufficiently differentiable 
conformal factor. The system (I2.16P has the form of Einstein equations with source 
terms. In general, such a system must be completed by additional equations derived 
from the Bianchi identity. The divergence freeness condition, V^T^i^ = 0, implies 
equations for the source functions. Some notable exceptions to this procedure are 
the Vaidya metric and the null dust [128j . In those cases, the divergence freeness 
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condition is satisfied without implying additional equations for the source functions 
so that certain functions can be prescribed freely, up to physically reasonable energy 
conditions, in terms of coordinates adapted to the null cones. It is interesting to note 
that these free functions are prescribed in a way that keeps the null cone structure 
invariant which is equivalent to the conformal structure. We show now that also in our 
case there are no additional equations required for the conformal factor f2 to satisfy 
the divergence freeness condition. We calculate 

+ -^{vM)iyxn)y^n. (2.17) 

Contracting the commutation relation, VaVj^ V^il — Va Vpfi = i?;^^p'^Vo-0, with 
g^P and exchanging derivatives we get, DV^ri — \7i,n\il — RJ'V ^^l. Using this, the 
definition of the Einstein tensor and the conformal source tensor as given in (|2.16p . we 
get the identity 

V^'T^u[n] = - ^ V^f} {G^,[g] - T^,[n]) . (2.18) 

We see that the divergence freeness condition, V''T^y[0] = 0, is satisfied identically by 
the conformally transformed Einstein equations for a non-vanishing f2 that is at least 
. Therefore, O can be regarded as a free function. We can write some favorable 
equation for its evolution consistent with the above calculation or prescribe it directly 
in terms of arbitrary spacetime coordinates as long as 17 7^ 0. 



2.3.1 A hyperbolic reduction 

We showed that the Einstein vacuum field equations for the physical metric g are 
equivalent to the conformally transformed Einstein equations (|2.16|1 for the rescaled 
metric g = fl^g with a non-vanishing, sufficiently differentiable conformal factor fl ^ 
written as a function of some yet unspecified coordinate system. However, this is not 
enough to work numerically with the system (|2.16p . We also would like to see that we 
can find solutions to (|2.16|1 . To show this we use the hyperbolic reduction technique 
[59]. 

The source terms in (|2.16|) involve second order derivatives of the conformal factor 
which translate to first order derivatives of the rescaled metric via Christoffel symbols, 

Vi,r2 — dfj^d^il — r^j^c^A ^- In a fully first order reduction of the Einstein equations, 
these Christoffel symbols can change the principal part. In the general wave gauge 
reduction we discuss below, however, the principal part does not change and the ar- 
gument for well-posedness from [69j can be directly used with a minor modification. 

Regarding the Ricci tensor R^''^ as a differential operator acting on the physical 
metric g, we can write the Einstein vacuum field equations in a local coordinate system 
{£^}p=o,i,2,3 as 

R'"'[9] - Ig^'dxdpr" + V(^f^) - -g^^r^ f^^f;, = 0, (2.19) 

where we have defined the contracted Christoffel symbols := g'^'^V'^^ = — Qgi^ 
and set V^f"" = gf'Pidpt'' + f^^^f^). The principal part of the operator R'"' is of 
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no known type. It was recognized by Choquet-Bruhat that one can always choose a 
wave gauge (historically referred to as harmonic gauge), at least locally, so that the 
contracted ChristofTel symbols vanish, F'^ = — Qgi^ — 0, and the Einstein vacuum 
equations l|2.19p . reduce to a quasi-linear system of wave equations. This reduction 
technique led to the first local existence result in general relativity [24] . 

The reduction based on the wave gauge was generalized to arbitrary coordinate 
systems by Friedrich with the introduction of gauge source functions |55j . In the 
general wave gauge, the coordinates are constructed as solutions to an initial value 
problem for the semi-linear system of wave equations Ogx'^ = —V^ = —F^^ with 
prescribed functions F^{x^'g) that can depend on the coordinates and the metric. 
These functions act as source functions for the coordinate gauge, hence the name 
gauge source functions. Note that the general wave gauge, in contrast to the wave 
gauge described above, is not a specific choice of coordinates but a way to prescribe 
general coordinates in an initial value formulation. 

The reduced system to (|2.16p is obtained by replacing the contracted Christoffel 
symbols with the gauge source functions F^^ . The resulting system becomes a quasi- 
linear system of wave equations for the metric components which can be written as 

G""'[c,] = r^'^[f]] + V^^C"^) - i(VAC•^)5^^ (2.20) 

where = F** — F^^ are called the constraint fields. This is a coupled system of 
quasi- linear wave equations for the unknown g^'' (|2.34p . We need to study the Cauchy 
problem for this system. I will just mention certain aspects that play a role in later 
considerations or that are different from the detailed discussion in ^69]. 

The Cauchy data on an initial hypersurface S = {x" — 0} consists of g^'^\s and 
do9^'^\s- Assume we are given on S a Riemannian metric haf3 and a symmetric tensor 
Kafi as a solution to the Einstein constraint equations and we have chosen gauge source 
functions F^''{x^). We further choose four functions on S that correspond to initial 
data for the lapse function a > and the three components of the shift vector . In 
the interior, these functions should be chosen such that 9o is timelike which implies 

— h^sP'^ > 0. Wc will later allow do to become null at the outer boundary (see 
the discussion leading to (12.32^ 1. We obtain the data g^'^\s via the decomposition 

g = g'^'d^.d, = -\dl + ^/3^9o5^ + ih'^' - ^) d^ds, j,S^l, 2, 3. (2.21) 

The data dog^'^\s is determined so that C^js = and Ka/s is the second fundamental 
form on S. Standard theorems guarantee that we can find a unique solution to the 
Cauchy problem for the reduced equations ()2.20|) that depends continuously on the 
initial data. Now we need to show that this solution is also a solution to (|2.16p . 
The solution spaces of (|2.20p and (|2.16p are equivalent if the constraint fields vanish, 
C = 0. We can derive a system of partial differential equations for the constraint 
fields by taking the divergence of our reduced system p.20p . The Bianchi identity, 
VpG^'' — 0, and (|2.18p imply the following subsidiary system for the constraint fields 

□C' + Ri^C - ^V^17 J^V^^'C"') - ^(VaC-^).?'"'^ = 0. (2.22) 

This is a homogeneous, semi-linear system of wave equations for . The Cauchy 
problem for this system has unique solutions. The initial data for the evolution equa- 
tions has been constructed such that the vanish on the initial surface, C^\s = 0. 
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From the evolution equations evaluated on S it follows that also the first derivatives 
of the constraint fields vanish on the initial surface, OqC^Is = 0. The uniqueness of 
solutions to the Cauchy problem for the subsidiary system (|2.22p then implies that the 
solution to the reduced system (|2.20p satisfies C' — away from the initial surface S. 
Therefore, we can use (|2.20p in a numerical calculation to find solutions to (|2.16p . 

The above argument does not change when we add homogeneous combinations 
of the constraint fields to the reduced system ()2.20p . These "constraint adjustment" 
terms can be of the form A'^{x,g,dg)C^ without changing the principal part of the 
reduced system and the uniqueness of the vanishing solution to the subsidiary system. 
Such terms can be useful for changing the behavior of the subsidiary system to damp 
constraint violations when the constraints are only satisfied up to an error. It is 
known that non-linearities in the subsidiary system can lead to catastrophic growth 
of constraint violations [66j . There are suggestions for constraint adjustment terms 
based on linearized studies |77| and numerical experiments with certain solutions [S]. 
A general procedure for their choice can hardly be expected. 

Note the appearance of a division by the conformal factor in the subsidiary system 
(j2.22p . With our current limited understanding of the subsidiary system away from 
the constraint surface, it can only be decided by numerical experiments whether this 
division causes difficulties for the propagation of constraint errors for a large class of 
dynamical solutions in regions where is small. 

We have seen that we can prescribe the representationjaf a conformal factor a priori 
in an initial value problem for the Einstein equations in A4, that is we can set D, = f{x) 
in (|2.20p where x'^ are coordinates in A4 and / is a positive function f{x) > 0, which is 
at least three times continuously differentiable. Such a prescription does not determine 
the conformal factor a priori as a function from the manifold to the positive real line. 
I shall point out how this prescription is to be understood to avoid confusion due to 
my sloppy but common notation. 

In an initial value problem we do not know the manifold a priori. The coordinates 
on the manifold are constructed during the solution process and they are determined 
by initial data as well as by the choice of the gauge source functions. The prescription 
of a function for a conformal factor determines only the representation of a conformal 
factor in terms of coordinates which are yet to be constructed during the solution of 
(|2.20p . Properties of the resulting conformal factor will depend on initial data and the 
choice of the gauge source functions for the coordinates, however, we can choose the 
coordinate representation of to be a convenient one for numerical calculations. 

The essential property of (|2.16p that is responsible for this feature is its con- 
formal invariance, in the sense that if {A4,g,n) is a solution to G'^jy[(?] = r^i/[r2], 
then (AijUj^g,Li!ft) with a sufficiently differentiable positive function w is a solution to 
G^^liv^g] = r^y[wl7]. The system (|2.16l) determines the conformal class of the metric 
g in contrast to the vacuum Einstein equations which determine the isometry class of 
g. In the spirit of Weyl, an arbitrary choice of a point-dependent unit of measurement, 
i.e. an arbitrary conformal gauge, is allowed. By prescribing a coordinate representa- 
tion for a conformal factor, we do not fix the unit of measurement, but its dependence 
on coordinates. The conformal factor is then constructed jointly with the metric, the 
manifold and the coordinates. 
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2.4 Including null infinity 

The freedom to prescribe a convenient representation for the conformal factor in terms 
of arbitrary coordinates is a feature valid for Q ^ where no invariant requirements 
are to be fulfilled by the conformal factor. However, our motivation to study conformal 
rescalings was to include J^"*" in the computational domain. To achieve this, we need 
to allow 17 = at the outer boundary of the future domain of dependence of an^jiiitial 
hyperboloidal hypersurface S while keeping fl > O^in the physical manifold A4. Of 
course, this can only be done, if at all, when {A4,g) admits a sufficiently regular 
conformal boundary. 

When is a prescribed function in the above sense, the requirement that {il = 0} 
shall correspond to J^^ implies certain conditions on the coordinates. The choice of 
the conformal gauge is then coupled to the choice of the coordinate gauge. In this 
section, we see how a suitable coupling can be achieved. 



2.4.1 The preferred conformal gauge at 

We assume that the Einstein vacuum field equations with vanishing cosmological con- 
stant are satisfied in a neighborhood of J^^. We restrict ourselves to the case where 
matter fields have spatially compact support, although a sufficiently strong fall-off be- 
havior of matter fields can also be regarded as compatible with our discussion. As 
a consequence of the Einstein vacuum field equations, ^+ is a shear- free null sur- 
face as discussed below (cf. |108) ). To find an appropriate conformal gauge for our 
calculations, we will use this property of 

Consider the transformation behavior of the Einstein tensor under conformal rescal- 
ings. Multiplying (|2.15p with fi^ and evaluating it along where fl = we see that 
g'^PV xil.'V pil\j!+ = 0. This shows together with dfl\j!+ ^ that is a null surface. 

We multiply (|2.15|) with fl and take its trace-free part along ^+ where i7 = to 

get 

1 



V^V.f} - -g^, m ) _=0, (2.23) 



The relation above is independent of the conformal gauge because we derived it from 
the transformation behavior of the Einstein tensor (|2.15p under some conformal rescal- 
ing. Another way to see this independence is to consider the transformation behavior 
of (|2.23p under a further conformal rescaling of the metric given by 

g' = Lo'^g, n'^ujVL, uj > Q on M. (2.24) 

Using dlJl) with V ^V + S{uj-'^duj), we derive the following transformation behavior 
for the second covariant derivative of the conformal factor along 

where we have set Vl\ = Va^I and ujx = The trace of the above relation reads 

D'n'\^+ = \ LjUVL + 4 VL^ijjy) . (2.25) 



We get 



V:,V' f]' - ig' , 



4 



= 0. 
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The relation (|2.23p can be interpreted as implying shear- freeness of J^+. To see 
this, we introduce in a neighborhood of ^+ a null vector field l'^ that satisfies on 
^+ the relation — V^il. We complete to a Newman-Penrose complex tetrad 
{l,k,m,fh) satisfying 

g{l, I) = g{k, k) = g(m, to) = g{m, rh) = 0, gil, k) = -gim, to) -1, 

gil, to) = to) = to) = m) = 0. (2.26) 

In [100] . Newman and Penrose introduce twelve complex functions called spin coef- 
ficients. We are interested in two of them defined by a :— to^to'^V^/i/ and p :— 
m^ifi^V fill,. As discussed in [lOOj . when l^^ is tangent to an afhnely parametrized null 
geodesic, a can be interpreted as the complex shear of the null geodesic congruence 
given by l'^ and the expansion of the congruence is characterized by p. Evaluating a 
and p at ^+ and using (|2.23p we see 

a\j;+ — —m'^m'^gfj^^Dn = 0, p\j'+ = — TO^TO'^g^i^DfJ = 

In our case, the null generators of are not necessarily geodesic, i.e. do not satisfy in 
general = on However, under a rescaling of l^^ given by (/')'' = Ol'^ with 

a positive function 9, the spin coefficient a transforms as a' = 9a, so the vanishing 
of a is invariant under a rescaling of l'^ which we can use to make geodesic. We 
therefore conclude that ^+ is a shear free surface. 

While the vanishing of a and the shear freeness of J^"*" is valid in any conformal 
gauge, the vanishing of the expansion of J^"*" characterized by □f2|j?+ depends on 
the conformal rescaling as seen from the transformation behavior (|2.25p . The relation 
(|2.25p implies that given fl whose expansion along ^+ does not vanish, we can always 
find a conformal rescaling (|2.24p such that □'r2'|j?+ = 0. To construct this conformal 
gauge in a given conformal extension, we need to solve an ordinary differential equation 
for the rescaling uj along the null generators of ^+ which reads 

This equation can always be solved in a given regular conformal extension and so, one 
can always construct the conformal gauge in which ^+ is manifestly expansion-free. 

We call the conformal gauge in which the expansion of ^+ vanishes, = 0, 

a preferred conformal gauge. This gauge is useful for analytic studies because of its 
special properties |1121 1129| . It is also used in numerical calculations in the character- 
istic approach [134j . A direct consequence of (|2.23p in a preferred conformal gauge is 
that = 0. By multiplying ((2T6l) with Q we see that in this gauge 

lim ^ g^'^VAflVflll = 0. 

So a useful property of this gauge is that each term in the conformally compactified 
Einstein equations (|2.16p attains a regular limit at J^^. 

Summarizing, a preferred conformal gauge in which ^+ is expansion free satisfies 
the following relations along 

an\jf+=Q ^ V^V^17U+=0 and lim ^g^'^VxQV [3^ = 0. (2.27) 
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2.4.2 Coupling the conformal and the coordinate gauge at 

In a preferred conformal gauge, we adapt the coordinates to J^^ so that ^+ is at a 
fixed spatial coordinate location. This secures that no resolution loss in the physical 
part of the conformal extension appears. 

Assume a conformal extension (Al,(7, 51) with a regular ^+ has been given in a 
preferred conformal gauge. As d£t\j:+ ^ 0, one can use the conformal factor as a 
coordinate in a neighborhood of which has the topology M x 5'^. We can then 
introduce another coordinate system {x^} — {t, r, a;"^} with respect to which SI = 1 — r 
in a neighborhood of Here, t is some time coordinate, r is a radial coordinate and 

are angular coordinates. This essentially corresponds to using the conformal factor 
as a coordinate in a neighborhood of J''^ which now corresponds to the coordinate 
surface {r = 1}. 

We use a specific representation of the conformal factor in terms of coordinates, 
ri = 1 — r, which couples the conformal and the coordinate gauge so that geometric 
properties of J^"*" translate into coordinate conditions. The relation SI = 1 — r satisfies 
(iri|j?+ 7^ 0. The requirement that SI is a preferred conformal factor such that the 
expansion of J^"*" vanishes manifestly in S7 implies 

= r''|r=i = 0. (2.28) 

By (|2.27p . we further have as a consequence of (|2.28p 

^ lim ig'^'^V.flV^ll = lim — = -drg"\r=i. (2.29) 

In this representation of the conformal factor, the source terms (|2.16p for the confor- 
mally transformed Einstein equations (j2.16p read 

T^.. = - Y^(r;. - 5M^r^) - (2.30) 

By our assumptions on the spacetime and our choice of coordinates, each term in T^y 
attains a regular limit along {r — 1}. 

We can specify dt further such that it is equal to = g^^Wiyfl along J^+. Then 
we get 

- -5''1r=i = = S^, (2.31) 

The relation (|2.3ip implies along = g""^ — g^^g^^ g^iu = gtt which is compatible 
with the fact that t is a parameter along a null surface. The vanishing of rj^^|r=i 
implies then in combination with (|2.3ip the following 10 relations along ^+ 

dtgtt = drgtt = dAgtt = 0, drgtr - ^dtgrr = 0, 
drgAt + dAgtr - dtgrA = 0, dAgst + dsgAt - dtgAB = 

The coordinate conditions we get are similar but more flexible than those in the char- 
acteristic approach [134] . This is because the hyperboloidal foliations allow for more 
gauge freedom. 
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2.4.3 Constructing preferred coordinates 

We have studied implications (|2.281 I2.29P of the specific prescription of a preferred 
conformal factor by = 1 — r in terms of a radial coordinate r. In a numerical cal- 
culation of a previously unknown spacetime, we would like to know how to construct 
preferred coordinates such that the coordinate surface {r = 1} has the desired geomet- 
ric properties of In general, it is not clear how to control geometric properties 
of given coordinate surfaces during time evolution because the spacetime which de- 
termines the geometry is constructed together with the coordinates, however, in the 
general wave gauge the condition (|2.28p can be controlled in a hyperboloidal initial 
value problem directly by an appropriate choice of the radial component of the 
gauge source functions. 

Assume hyperboloidal initial data is given on a spacelike hypersurface S with a two 
dimensional boundary E such that its evolution admits a regular in accordance 
with the theorems of ^3, i54j . We can always do a coordinate transformation on S such 
that the conformal factor has the form f^ls^l — rina neighborhood of S. 

In the hyperboloidal initial value problem for the conformally compactified Ein- 
stein equations using the general wave gauge reduction (|2.20p . the coordinates are 
constructed during time evolution implicitly as solutions to the wave equations with 
source terms □gx'^ — — = — F^. As seen by this relation and = 1 — r, the radial 
gauge source function F*" at {r — 1} determines directly the expansion of We 
can make sure that (I2.28P is satisfied by simply setting = initially and during 

time evolution. The initial data for lapse and shift should be chosen according to 



By our choice of initial data and gauge source functions we have a preferred conformal 
gauge and therefore also l|2.29p . Note that while it is not clear whether (|2.3ip can be 
satisfied by the choice of gauge source functions, (|2.29p is sufficient for the conformal 
source terms given by (I2.30p to attain regular limits at ^+ that can in principal be 
calculated numerically. A possibility to deal with these formally singular terms is 
discussed in section [231 

We note again that, even in a given manifold, a prescription such as J7 = 1 — r will 
result in different functions on the manifold depending on the choice of the coordinate r. 
It is a remarkable feature of the conformally compactified Einstein equations that they 
allow us to prescribe the conformal and the coordinate gauge in a way that geometric 
properties of are respected by the coordinates in an initial value formulation and 
the representation of a preferred conformal factor in terms of coordinates is known a 
priori. 

2.4.4 Choice of a conformal factor in the interior 

The notion of an isolated system implies the existence of an interior and a far-field 
region. The suggested method to include null infinity in the numerical domain is 
tailored for the treatment of the asymptotic region. The numerical evolution scheme 
in the interior does not need to be changed in_this approach as we are free to prescribe 
the representation of a conformal factor in M.. 



([23T|) . By (f23T|) we set 





(2.32) 
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The interior domain on a spacelike surface will be confined to the interior of a ball 
of radius which includes the spatially compact support of matter sources. Assume 
that a radial coordinate r has its origin in the center of this ball and choose a large 
ra > Ti. The domain < r < can be regarded as the transition domain between 
the interior and the asymptotic region. The conformal factor in the interior can be set 
to unity, U, — 1 for r < fi, so that the conformal source terms in (I2.16|) vanish. In the 
asymptotic region r > we can choose il = i? — r with any large constant R > ra 
without changing the analysis in the previous sections made for a neighborhood of J^^ 
which is now given by {r = R}. On < r < 7'a, any transition between the different 
domains can be used that is at least C^. A smooth transition is given for example by 



There are of course many other possibilities. The important point is that we are allowed 
to prescribe a non-vanishing, sufficiently differentiable function for the representation 
of a conformal factor to solve the Einstein equations in the interior which allows us in 
principle to attach a compactified asymptotic region to standard numerical rclativistic 
calculations. The width of the transition region or the form of the transition function 
should be decided upon by empirical studies. 

2.5 Numerical tests in spherical symmetry 

A question that we did not discuss so far is how to set up a numerical code such that 
the formally singular source terms T^i, as given in (|2.30p can be dealt with. In this 
section, we study a simple possibility and discuss numerical test results on the example 
of the extended Schwarzschild spacetime as given in (|2.13p . 

The system (|2.20p is not suitable for analytic studies on questions of smoothness 
properties of ^+ , as the equations are formally singular at the set whose properties we 
are interested in. In this section, however, the main interest does not lie in questions 
of differentiability or existence but in the numerical construction of spacetimes which 
admit a sufficiently smooth conformal boundary. To achieve this in the case of spherical 
symmetry, we numerically solve the conformally transformed Einstein equations in 
the general wave gauge discussed in previous sections for the extended Schwarzschild 
spacetime. The example we study is clearly very special. 

The general wave gauge has been used in various numerical studies [71 [51 [7^ [Ml 11061 
HHinZlinSKIlllinilllSKnS]- There are many aspects to numerical evolutions that 
a successful calculation needs to deal with. The numerical results presented below are 
not optimal with respect to the numerical grid boundary treatment or the constraint 
propagation. The aim is only to see whether the suggested method is robust enough 
to deliver results with a simple choice of evolution variables and a straightforward 
treatment of the outer grid boundary. 



n = e ^^"^^ + {R-r)[l-e (^^^ 
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2.5.1 The evolution system 

The evolution system can be written as 

RAa] - ^ipC,} + = T^A^] ~ i g^.Tp], (2.33) 

with — i^p. Here, A^j^^ are constraint adjustment terms of the form 

A^u = CxA^^{x,g,dg). 

These terms vanish when the constraints are satisfied. They change the propagation 
properties of constraint errors without effecting the principal part of the system. In 
terms of metric components we get for (|2.33p 

lo^'dxdpg^, - V(^F,) - rl^T^.g^'gna - 2r^,5^''5,(^r:^)^ - A^, = 

= (2r;. + F^g,.) - jrhr'^^'"'- ^^-^^^ 

For the right hand side, = 1 — r has been set explicitly and the contracted Christoffel 
symbol has been replaced by the gauge source function F"^. This replacement 
modifies the subsidiary system (|2.22p by a single term without changing the results 
we discussed in 12.3.11 

The calculations are done for a spherically symmetric metric in adapted coordinates 
(t, r, 'd, if). The line element can be written as 

g — gtt dt^ + 2gtr dt dr + grr dr^ + g^-a da^ , (2.35) 

where da^ = di?^ + sin^-d d(p-^. All metric components are functions of {t,r) only. 
Alternatively one can write the metric in terms of the variables of a 3+1 decomposition 
as 

g = {-a^ + h^p^)dt^ + 2h^pdtdr + h^dr^ +b^da^, 

where the lapse a, the shift /3, and the spatial metric functions h and b arc functions 
of the coordinates {t, r) only. They are given in terms of the metric components by 

1 gt n 9tr , . 

« = . jf = \ gtt, P = , /l = V5rr- 

V-ff V 9rr grr 

The evolution variables can be chosen such that the resulting system of partial 
differential equations is first order in time and second order in space. This can be 
done in various ways with different stability properties. A simple choice is given by 
g^i, and the auxiliary variables tt^i^ = n^Oxg/^n, where n'^ is the unit normal to the 
surfaces t = const, given by n = {dt — Pdr). We write the system of evolution 
equations for the variables (g^i/jTr^i/) as 

dtgp,v = P drgp^u + aTTf_,„, 

StTT^i/ = fidriT^u + — dlgf_,u + H^^i,{g,drg,TT,r), 

grr 

where the H^^^ are lower order terms. The initial data at {t = 0} is read from the 
explicit solution (|2.13p and no boundary data is needed. 
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The numerical code is based on finite differencing for the spatial derivatives and 
method of lines for time integration. We discretize the computational domain by 
introducing the homogeneous grid rj = (tq + jSr). The radial coordinate location 
of the spacelike inner boundary inside the event horizon is denoted by tq, the grid 
spacing is denoted by Sr and j — 1,2, . . . ,n where n + 1 is the total number of nominal 
grid points. We define difference operators D± by their action on a grid function 

via D^(l)j = — (pj)/{6r) and D^cfij — {(pj — 0j_i)/((5r). We also define 

Dq = + Z3_)/2. The spatial derivatives are replaced by fourth order accurate 
discrete derivatives 



drcb^Doil- ^ D+d\ 0„ ^ D+D- (1 - ^ D+D^ ) cj),. (2.36) 



The calculation of derivatives at the boundaries is described in l2.5.3l To each evolution 
equation for a variable (j) we add Kreiss-Oliger type dissipation terms given by 

where e is a small dissipation coefficient. No dissipation is applied in a neighborhood of 
the boundaries. The time integration is done by a fourth order accurate Runge-Kutta 
scheme. 



2.5.2 Choice of the gauge source functions 

For the evolution of Einstein equations starting from some initial surface, we prescribe 
the covariant components of the gauge source functions, The form of the metric 
(|2.35p and the Einstein equations (j2.33p imply for the angular components 

F^^- cot z9, = 0. 

The gauge source functions Ft and Fr are free. In spherical symmetry they can be 
prescribed in a way that distinguishes the radial coordinate r. For any metric adapted 
to spherical symmetry, one can find a radial coordinate such that g^^{t,r) — or 
b{t,r) = r. Such coordinates are called ureal as in these coordinates any r ~ const. 
2-surface on a time slice has the area Airr'^. 

In the ADM formalism, the condition that r is an areal coordinate is called the " area 
locking" condition. It determines the shift by an algebraic relation ^5S]. In the general 
wave gauge, the condition that r is an areal coordinate implies a relation between 
the gauge source functions Ft and Fr. This relation can be derived from the Einstein 
equations by setting g^i){t,r) = or b{t,r) = r, into the angular component of the 
equations. Without compactification {T^i, = in (|2.33p ) and constraint adjustment 
{Afii' = 0), we get the following relation between the gauge source functions Ft and Fr 

9tt - 9tr + 9tt9rr + r {Fr Qtt " Ft Qtr) = 0. 

This relation can be used to determine for example Ft{gtt, 9trT grr, Fr,r) for a given 
Fr- In terms of variables of the 3+1 decomposition we get 
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for non- vanishing shift. For the conformally transformed Einstein equations one gets 
the relation 

F, = El9ll + J_ ((_ 2 ^ g^^g^^^^ ^ g^^^ _ ^ (2.38) 
9tr gtrr gtAl 

Unfortunately, the areal condition in this form is not suitable for our calculations in 
general, because the shift (or equivalently gtr) vanishes in the interior. This is because 
we require the shift to point into the computational domain at both boundaries. For 
the interior calculation around the excision region of the black hole, the shift is required 
to be positive so that the inner boundary is a spacelike surface inside the event horizon. 
At the outer boundary, the shift becomes negative since is an ingoing null surface. 
The change of sign in the shift seems to be a general feature of hyperboloidal foliations 
that can be used for excision. 

For the numerical tests, we do not use (|2.38p . We read the gauge source functions 
from the explicit solution (I2.13p . In spherical symmetry, the condition that = 

corresponds to Ft\y+ = 0. A simple calculation shows that this condition is satisfied 
in our case. 

2.5.3 Numerical treatment of grid boundaries 

We introduce ghost points for the discussion of grid boundaries. To build fourth order 
finite differences via (|2.36|) at the boundaries we need two ghost points. Our method 
for the numerical treatment of boundaries for a first order in time second order in 
space system relies on [28] . Below, we give prescriptions for the grid functions at the 
inner ghost points j = —2, —1 and at the outer ghost points j = n + \,n + 2. For the 
simulations, we choose m = 1 and set the boundaries at 

2 2 1 

''0=0 ' rn = l , 

in n 

which implies a grid spacing of dr — (r„ — ro)/n = l/(3n) — The event horizon 

is at 2/3 and ^+ is at r = 1. For large n, the distance of the inner boundary to the 
horizon is about 6 (Sr) and the distance of the outer boundary to is about 3 {6r). 
Note that the coordinate distance decreases with increasing resolution. 

Inner Boundary 

The inner boundary is a spacelike surface inside the event horizon of the black hole as 
depicted in the Penrose diagram Fig. 12.71 or in the diagram of the causal structure on 
the grid Fig. 12.51 No boundary data is needed on a spacelike hypersurface so that we 
can excise the singularity from the computational domain. The computation outside 
the event horizon should not be disturbed by the excision inside the black hole. 

The numerical outflow boundary conditions consist of a fifth order extrapolation 
for g and fourth order extrapolation for tt [55]. We set 

{SrfDlg.i=0 => g_i = 5go - lOgi + lO.ga - 5^3 + 34, 

iSr)^Dlg^2 = => g_2 = 5g_i - IO.90 + IO51 - 5^2 + 53, 

((5r)^_D^}L7r_i =0 =^ 7r_i = 4 tto — 6 tti + 4 7r2 — tts, 

{Sr)^D'lTT_2 ^ 7r_2 = 4 7r_i - Gtto + 4 7ri - 7r2. 
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Fig. 2.6: Stencil for extrapolation at the outer boundary. 



Outer Boundary 

There are no characteristics entering the computational domain from the outer bound- 
ary therefore no analytic boundary data are needed. A numerical treatment, however, 
is necessary to deal with the appearance of divisions by factors of f2 on the right hand 
side (I2.30p . Two possibilities are: 

• Calculate the limit for T^j^ at by numerical techniques (for example by using 
a discrete version of the I'Hospital rule) . 

• Extrapolate the solution (g^jy, tt^;/) at the outer boundary. 

I chose the second option, not because it is a clean treatment, but because it is a 
simple method to implement. In the simulations presented below, .y^ is about three 
grid points away from the outer grid boundary (Fig. 12.6)) . The ghost points are filled 
by the outflow conditions as for the inner boundary 

{SrfDtgn+i^O .gn+i = 5.g„ - 10.g„_i + 10.g„_2 - 5g„_3 + .g„_4, 

^ .gn+2 = 5 .g„+i - 10 .g„ + 10.g„_i - 5 .g„_2 + ffn-s, 

((5r)*I?lTr„+i =0 TTn+i = 4 7r„ - 6 7r„_i + 4 7r„_2 - 7r„_3, 

((5r)''£)i7r„+2 = => 7r„+2 = 4 7r„+i - 6 7r„ + 47r„_i - 7r„_2- 

Note that in this method, the outer grid points move along a timelike curve and is 
not on a grid point. Clearly, this method is not optimal. Extrapolation near .y^ is not 
very accurate. An accurate calculation in a neighborhood of is however important 
for a general calculation of gravitational radiation that depends via the rescaled Weyl 
tensor on third derivatives of the metric. To establish the idea presented in 12.31 and 
12. 4[ a more sophisticated numerical treatment of the formally singular terms will be 
required. As mentioned before, our aim in this section is not to construct a general 
numerical method, but to test the suggested idea on a special example using simple 
techniques. 



We test our method on the example of the extended Schwarzschild spacetime in a 
CMC foliation presented in 12.2.21 Malec and Murchadha discuss in the asymp- 
totic behavior of the embedded CMC slices depending on the parameters K and C. 
The surfaces K > reach future null infinity. Their behavior in the interior depends on 
the choice of C. Below the critical value C = SKm^/S the surfaces pass the Schwarz- 
schild throat below the bifurcation sphere. Fig. 12.71 shows foliations of the extended 
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Schwarzschild spacetime with the same value of C = 2, m = 1 and three different val- 
ues oi K = {1,0.3,0.07}. In the case K = 1 the surfaces go through the past horizon 

into the white hole. Then the shift vector given in (j2.14|l by f] — — '^^^^ ^2^'^^ ct is 
negative over the computational domain and excision can not be expected to work. 
In this case, no numerically stable evolution was possible. Note that in a black hole 
spacetime formed by gravitational collapse, no bifurcation sphere and no past horizon 
exists. 

Below we discuss numerical test results for the cases K = 0.3 and K — 0.07 where 
the condition C > 8Km^/3 is satisfied with C — 2 and m — 1. As seen in Fig. l2.7l these 
surfaces go through the future horizon. The dissipation coefficient has been chosen to 
be e = 0.05 and no constraint adjustment terms have been added. The solutions have 
been calculated on the domain r G [2/3 — 2/n, 1 — l/n] where rt + 1 is the number 
of grid points. We use 501, 1001 and 2001 grid points on a homogeneous grid. In 
terms of the Schwarzschild radius f, the lowest resolution correspond to the domain 
r G [1.96to,499to] while the highest resolution corresponds to f e [1.99m, 1999m]. 

Fig. l2.8l shows the convergence factors c(t) in the _L2-norm for the metric component 
gtt as a function of time. The convergence factor c{t) for a grid function (j) is calculated 

by 

C{t) ^ l0g2 



where cj/^'^'^ is the numerical solution in medium resolution, (/I'^'f^ in high resolution 
and (j)'^'^ is the explicit solution. The L2 norm of a real function (f) £ i^(iSt, M) over the 
one dimensional computational domain St reads 

^dr]\ (2.39) 
St / 

The discretized version of the L2-norm reads 

(2.40) 

To compare different resolutions for the convergence tests, the sum over the grid points 
in the L2 norm is sampled according to the resolution. 

We observe in Fig. 12.81 for K = 0.3 that the convergence factor drops to roughly 3 
after about 40to. For K = 0.07 the decrease is faster. A global convergence factor of 3 
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might be expected due to our numerical outer boundary treatment. After about 80to 
the codes crash. To localize the problem of the crash, we plot in Fig. 12.91 and Fig. 12.101 
the radial dependence of the error for gu at times t = {IOto, 20to, 30to, 40to}. The 
aim is to see on which part of the computational domain the errors become large and 
where convergence is lost. 

For representing the errors of different resolutions in the same plot, we scale them 
depending on a convergence factor c. For example when the convergence factor is 
expected to be 4, the error for the solution in medium resolution is scaled by 2'' = 16 
and the error for the solution in high resolution is scaled by (2^)'* = 256. When the 
curves are on top of each other, we can assume that the overall convergence factor is 
close to the chosen one with which the errors have been scaled. This visual test gives 
us a possibility to localize deviations from convergence along the grid. 

In Fig. 12.101 we see that already at t — 20m the convergence in the interior is dis- 
turbed. At t — 30m the convergence factor has dropped to about 3.1 while the lowest 
resolution does not align with the medium and the high resolution. This deviation of 
the low resolution can be seen more clearly at 40to. For K = 0.07 the largest errors 
appear in the interior where the code eventually crashes. 

In Fig. 12.91 with K = 0.3 we see that the decrease in the convergence factor is 
slower in accordance with Fig. l2.8l but this time the largest errors are due to the outer 
boundary treatment. One observes deviations from convergence in the interior too but 
the error is dominated by the region close to = {r = 1} and the code crashes in 
this region. 

For small K the numerical errors produced in the interior propagate outwards 
slower compared to large K as can also be seen in the plots of causal structures 
Fig. 12.51 Therefore the numerical errors in a neighborhood of J^"*" stay small for small 
K while the errors in the interior grow until the code crashes in contrast to large 
K which crashes close to J^^. We observe that different choices of K have a strong 
effect on numerical errors. A related effect of the foliation has been studied in [S2] 
where it was shown for a certain system that the choice of the time slicing influences 
propagation properties of constraint errors. 

Summarizing, we can say that although our numerical setup does not allow us to do 
long time evolutions of the extended Schwarzschild spacetime, a piece of null infinity 
can be calculated with the method suggested in 12.31 and 12.41 even with simple choices 
of variables and numerical boundary treatment. 



2.5 Numerical tests in spherical symmetry 



48 



t=10m, c=4.1 



t=20m, c=4 





8e-09 r 




6e-09 - 






O) 
_g 


4e-09 - 


rror 


2e-09 - 


LU 






Oe+00 - 




-2e-09 - 





t=30m, c=3.7 



t=40m, c=3.3 



1e-05 
8e-06 
6e-06 
4e-06 
2e-06 
Oe+00 





low ' — 




2= med ------ J 




2^'= high 

J 


















0.7 


0.8 0.9 




Fig. 2.9: Rescaled error plots for K = 0.3 
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Fig. 2.10: Rescaled error plots for K = 0.07 
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2.6 Discussion 

The main content of this chapter has been the development of a numerical method 
that allows us to construct solutions to the Einstein equations including null infinity in 
the computational domain based on a general wave gauge with ^+-fixing coordinates. 
The method employs a suitable coupling of the conformal and the coordinate gauge 
to establish expansion freeness of J^"*" via an appropriate choice of the gauge source 
functions for the coordinates. Under our assumptions on initial data and our gauge 
conditions, each formally singular term arising from the conformal compactification 
attains a regular limit at which needs to be calculated by numerical techniques. 

I presented a numerical test of this method in spherical symmetry in which the 
outer boundary treatment is based on extrapolation near J^~^ . The code is afflicted by 
an unbounded growth of numerical errors. One could see, however, that the treatment 
of the outer boundary does not lead to immediate problems so that a piece of null 
infinity can be calculated. It is an outstanding question whether one can calculate 
the formally singular terms arising from the conformal compactification numerically 
in a stable manner for dynamical spacetimes. To achieve this one might need to use 
intrinsic properties of the characteristic surface to a larger extent. 

As a next step, one should try the method described in 12.31 and 12.41 in a setting 
including radiation using a more sophisticated choice of variables and numerical bound- 
ary treatment. Techniques developed in various groups over many years working with 
the general wave gauge reduction might be useful in that context ^ J35l I106[ I114[ I122[ 
I131j . One should also study the freedom in prescribing gauge source functions further 
to derive useful gauge conditions which lead to a convenient representation of the so- 
lution metric along .y^ . A numerical finite differencing scheme might be adapted to 
certain expected properties of the solution metric as in the characteristic approach. 

As has topology M x 5'^, a general numerical calculation in which is fixed 
to a spatial coordinate location which corresponds to the outer grid boundary should 
be able to handle spherical grid topology. To my knowledge, there are currently two 
approaches in numerical relativity that avoid coordinate singularities on a sphere and 
use the general wave gauge |106[ I122j . Both of these methods seem promising for 
trying the suggested idea in a general setting. 

In our studies we used a general wave gauge reduction. However, the idea to solve 
a hyperboloidal initial value problem for conformally compactified Einstein equations 
with a prescribed representation of the conformal factor is independent of the reduc- 
tion. The open question for other reductions of Einstein equations is whether the 
conditions (|2.281 12.29P can be satisfied during time evolution in a well-posed Cauchy 
problem. This question needs to be studied for each reduction separately. 

The maximal development of hyperboloidal initial data does not yield the global 
spacetime. We do not get access to spatial infinity in a hyperboloidal initial value 
problem and in J^+-fixing coordinates, also timelike infinity can not be reached. While 
one can argue that these points are not of physical interest and possibly irrelevant 
for the comparison of observational data with numerical calculations, it would be 
desirable to have access to the global structure for various reasons. In the hyperboloidal 
approach it is not clear whether the cut of the initial hyperboloidal surface at null 
infinity is close to timelike infinity or to spatial infinity. One would also like to be 
able to relate asymptotic quantities as mass or momentum defined at null infinity to 
corresponding quantities at spatial infinity. 



Chapter 3 

Spatial Infinity 



In this chapter we want to solve numericaUy a Cauchy problem for the Einstein equa- 
tions starting from a Cauchy hypersurface including spatial and null infinity in the 
numerical domain. This gives us, in principal, access to the global spacetime solution. 
We will see that in spherical symmetry the maximal development of Schwarzschild- 
Kruskal initial data given on a Cauchy hypersurface can be calculated. For more 
general situations including gravitational radiation we will focus our interest on the 
calculation of the detailed structure of gravitational fields in a neighborhood of spatial 
infinity including a piece of null infinity. 

We make certain assumptions on initial data that allow us to discuss the main fea- 
tures of our problem while simplifying the calculations involved. We consider asymp- 
totically fiat, time symmetric, vacuum initial data, i.e. we are given (5, ha/3), where S 
is a three dimensional manifold with an asymptotically fiat end and hap is a positive 
definite, asymptotically fiat, Riemannian metric on S. The data satisfies the vacuum 
Einstein constraint equations with Kap = 0. We mean by asymptotic fiatness for 
(5, hap) that the complement of a compact set in S is diffeomorphic to the comple- 
ment of a closed ball in and in the chart {x"} given by this diffeomorphism the 
following fall-off conditions are required 



with e > 0, fc > 2 and m is the ADM- mass of the initial data set. The fall-off conditions 
above are written with respect to a coordinate system. To discuss such conditions 
without relying on specific coordinate systems one can use conformal techniques along 
the lines of Penrose's considerations. In [71] Geroch uses conformal compactification 
techniques to represent spatial infinity as a single point i so that asymptotic properties 
of fields on S may be treated as local geometric properties at i. The construction is 
similar to the conformal mapping of R"^ onto the 3-sphere by adding a single point at 
infinity. We require that there exists a manifold S ^ S U {i} and a conformal factor 
with (/) > on S satisfying 

• h — (jy^h is & smooth metric on 5, 

• = 0, Da4) — 0, DaDpcj) — 2 hap at i. 




00, 



a, (3 — 1, 2, 3, 
(3.1) 
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where D is the covariant derivative operator defined by h. The described notion of 
asymptotic flatness is a property of the conformal structure in the sense that the 
conditions are invariant under conformal rescalings h — > uj'^h, (p —> uxj) with positive 
smooth functions lo satisfying uj{i) — 1. In the examples that we will study, h will not 
just be smooth but analytic at i. 

To calculate numerically an entire, asymptotically flat spacetime including spatial, 
null and timelike infinity, we might wish to solve Friedrich's conformally regular field 
equations presented in [53j . The reason why we can not do this is that the initial 
data for the conformal field equations blow up at spatial infinity when represented as 
a point i for non-vanishing ADM-mass. Specifically, denote by r(p) the distance of the 
point p €z S from the point i in terms of the metric haf3 on the Cauchy surface S. It 
turns out as a consequence of the constraint equations implied by the conformal field 
equations on S that the rescaled Weyl tensor blows up like [60] 

W)ki = O as r ^ with m^O. (3.2) 

A regular, finite, initial value problem near spatial infinity could be formulated by 
Friedrich based on an extended system of conformal field equations and a certain con- 
formal gauge near spatial infinity, called the conformal Gauss gauge. This gauge allows 
a representation of spatial infinity as a cylinder denoted by X. In this representation, 
the point i is blown up to a sphere (Fig. 13. ip . The blow-up procedure is such that 
a suitable rescaling of fields on S results in smooth conformal data up to and beyond 
1°. The development of this data can be studied further by the evolution equations. 
A basic difficulty in this new representation is the degeneracy of the equations at the 
sets where null infinity meets spatial infinity. 

The construction of the cylinder at spatial infinity is a delicate procedure which 
relies on a deliberate choice of the coordinate, the conformal and the frame gauge 
in accordance with a suitable notion of asymptotic flatness. For further details of 
this construction that we do not discuss in this chapter, the reader is referred to the 
original paper [ST] and the review articles |SH [M] . Further applications of the cylinder 
at spatial infinity in analytic work can be found for example in [1391 11421 1143| . 




Fig. 3.1: Point compactification and the cylinder at spatial infinity. 

We are concerned in this chapter with a numerical implementation of the extended 
system of conformal field equations with the cylinder at spatial infinity. First we 
discuss the numerical construction of a conformal Gauss gauge in the simple examples 
of Schwarzschild and Kerr spacetimes. After discussing some aspects of the reduced 
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general conformal field equations, we integrate the equations in spherical symmetry 
for the Schwarzschild-Kruskal spacetime. To include radiation into our discussion 
in a simple setting in which the problem (|3.2p does not appear, we calculate initial 
data with vanishing ADM-mass but a non- vanishing radiation field along J^+. The 
numerical development of this data is discussed in a neighborhood of spatial infinity 
represented both as a point and as a cylinder. The numerical implementation of the 
equations with the cylinder is such that the code can also be used to study physically 
more interesting spacetimes with non-vanishing ADM-mass. Such studies are left for 
future work. 

In the following, we will use a frame formalism. Latin letters are used for frame 
indices with i, j,k, . . . = 0, 1, 2, 3 and a,b,c, . . . = 1,2,3. Greek letters are used for 
coordinate indices with /x, i/. A, . . . = 0, 1, 2, 3 and a, /3, 7, . . . = 1, 2, 3. 

3.1 The conformal Gauss gauge 

The main reference for the properties of the conformal Gauss gauge that we discuss 
below is [63j. This gauge is based on conformal geodesies. 

3.1.1 Conformal geodesies 

Null geodesic congruences, when they are smooth, provide a valuable tool to study 
the asymptotic and causal structure of spacetimes. As null geodesies are invariants of 
the conformal structure, one might presume that spacelike or timelike curves that are 
conformal invariants might also be useful in such studies. In general, geodesies with 
respect to a metric are not geodesic with respect to a conformally rescaled metric. 
Conformal geodesies, however, are conformally invariant in the sense that, as point 
sets, they are independent of the metric chosen in the conformal class. They are 
autoparallel curves with respect to a Weyl connection defined in 11.11 

A solution to the conformal geodesic equations given below does not only provide 
a spacetime curve, but along the curve also a Weyl connection, a conformal factor and 
a frame which is orthonormal for a metric in the conformal class. While the equations 
are independent of coordinates, we will have to write them in some coordinate system, 
as we will study them in numerical applications. 

Given a metric g, the equations that define a conformal geodesic x'^(t) are written 
for its tangent vector x'^{t) and a covector /^(t) as 

{^ixr + s{f)^^^x^ip = 0, 

(V±/),, - i/A^(/W = (3.3) 

where S{f) has been given in (|1.2p and = ^ R^i, — Rg^iv is called the Schouten 
tensor. It is related to the Schouten tensor L^j, of a Weyl connection V = V -I- S{f) 
by i^i/ = L^^ - V^/i. -f f^fv - ^gt^ug^Pfx/p- By setting = in we see that 
every physical geodesic of a vacuum spacetime is also a conformal geodesic. 

We can write the equations (|3.3p using the Weyl connection V that is defined by 
the 1-form / along the curve x^(t) by V = V -|- S{f). The equation for x^^ becomes 
Vii^ = 0, so that the curve a;^(r) is autoparallel with respect to the Weyl connection 
V. The equation for / becomes i^'i^jy — 0. We will also use a frame that is V-parallelly 
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transported along the conformal geodesies. The equation for the components of the 
frame vector fields ej,''(T) reads '^x^k ~ ^- ^'^^^ take x to be the timelike frame 
vector field eg. 

We see by that the metric g = fl^g is V-parallelly transported along a given 
curve x(t) in M if the function > satisfies on x{t) the equation 

yin = nf^xr (3.4) 

In a spacetime solution to the vacuum Einstein equations one can derive for the confor- 
mal factor il(r) satisfying (|3.4|) the equation V|il — (see [S^ or [63]). This equation 
can be solved explicitly so that the conformal factor ^(t) is known a priori in terms 
of initial data. This is a remarkable property of conformal geodesies which will play 
an important role in this chapter. For given initial data on 5 = {t = 0} 

n{o)^n,, x^(o) = xj:, i'^(o) = ii:, /^(o) = (/,)^, 

the conformal factor 0(r) along the conformal geodesic can be written explicitly as 

^{t) = n.(^l + T {Uxn\r=o + (5M.i'^i'^)|r=o (s^Va/p) k=o) . 

It is useful to distinguish between the conformal compactification of the induced 
metric h on the initial Cauchy hypersurface S* = {r = 0} and the conformal compact- 
ification of the spacetime metric g. We write for the compactification of the induced 
spatial metric h — (p'^h, with (p satisfying the relations given in the beginning of the 
chapter. The choice of fJ* can be made such that spatial infinity is represented either 
by a point or by a cylinder. We introduce a free function k by setting = -. When 
we choose K\i = 1, dK\i = 0, then spatial infinity i corresponds to the point z° of the 
spacetime. Choosing /t|i = 0, duli ^ 0, results in the representation of spatial infinity 
as a cylinder X (Fig. 13. ip . We will make use of this freedom in later sections of this 
chapter. 

We have h — VLih — ^rh — -Kli. We choose initial data for the conformal geodesies 
such that (/*)^a;* = 0. On the initial hypersurface S we set [f^.)^ — By 
requiring .g(i,i)|r=o = ~lj we have g{x,x)r={) = Then the conformal factor 

becomes 

^^r)^'t(l-r^^\ with c.^ ^'^ =■ (3.5) 

Another important role of the free function k is seen by this formula. We can control 
the value Tjr± = ±^ of the time coordinate at which the conformal geodesic cuts 
= {f](rjr±) = 0} by the choice of k. 
Beside the conformal factor, the field dk '■— ^ ftJ-^k '-^^ determined explic- 

itly. It turns out that 

dk = mT),dam = f-2r^,i(e,)/a^0) , (3.6) 

where (e,)^'' :— ej^{0). For a derivation of these results, see [SH]. The knowledge 
of the conformal factor by p.5p will be very useful in our numerical studies in later 
sections. 
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3.1.2 Construction of the conformal Gauss gauge 

To construct a conformal Gauss gauge one uses conformal geodesies in a similar way 
as one uses metric geodesies to construct the Gauss gauge. One specifies a congruence 
of timelike vectors, a conformal factor and a 1-form on an initial hypersurface. The 
timelike conformal geodesies starting from this surface provide the conformal Gauss 
coordinates. Spatial coordinates are dragged along. 

Assume a solution {M , g) to the vacuum Einstein equations has been given. We 
construct a conformal Gauss gauge with timelike conformal geodesies and an orthonor- 
mal frame along them, (x'^(t), /^(t), eJ^{T) ), as follows: 

1. We find a conformal extension {A4,g,9) of the solution {Ai,g). In general, we 
will need different coordinates in different asymptotic regions which are used for 
the calculation of conformal geodesies via (|3.7p . 

2. On a spacelike slice with a spatial metric h, we introduce compactifying co- 
ordinates and rescale h with a suitable conformal factor 0, such that in these 
coordinates we have h = (jx^h. The metric h is used in the calculation of the 
conformal factor p.Sp . 

3. We set initial data (x* , iC, (e*)^) = (x^, i^, e^)|s according to 

{f*)t^ = (t> ^df,(j), (/^i'')|s = 0, _L S", g{x,i)\s = —gix,i)\s = -1, 

5(ea,i)|s = 0, g{ea,eb)\s = Sab, where a, 6=1,2,3. 

The timelike frame vector is given by x itself. The frame is not unique, the 
freedom corresponds to the freedom of spatial rotations. One also needs to 
choose the free function k that determines the representation of spatial infinity 
and the value of the time coordinate on in the conformal Gauss gauge. 

4. We solve the following system of ordinary differential equations 

{drxT - -f /^i^i" - 2i,Ux'') x'^ + {gxpx^x") 

idreaY = -V^^x^eP - {U:) - {li'^X + {dxpe^i") (3.7) 

The one-form / for which we solve the equations is related to the Weyl connection 
V = V + S{f) by V = V + S{!) where V = V + S{B-^dQ), so / = / + Q-^dQ. 
The initial data for / is chosen accordingly as /* = {(f) / 0)'^ d{(t) / 0) . 

In a numerical calculation, we check the quality of the solution using (|3.5p and (|3.6p . 

We introduced various conformal metrics to do the calculation. The physical metric 
is denoted by 5, it induces a spatial metric h on the initial hypersurface S which we 
compactify on 5* = 5 U {i} by rescaling h = 4>^h. As the conformal geodesies cover 
asymptotic regions where the physical metric g becomes singular, we also introduced 
a compactified spacetime metric g = O^g to calculate the right hand side of (|3.7p . 
Another conformal spacetime metric is the one that is V-parallelly transported along 
the conformal geodesies with the conformal factor satisfying p.4p . It is acquired by 
g = ri^g with Vl given explicitly by (|3.5p . 
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3.1.3 Numerical experiments 

In [63j Friedrich analytically constructs conformal Gauss coordinates in the Schwarz- 
schild-Kruskal spacetime and shows that they cover the conformal extension in a 
smooth way. In the following, we will construct such coordinates numerically. Go- 
ing beyond the analytical studies, we will also see in our numerical studies that one 
can construct conformal Gauss coordinates on the Kerr spacetime. 

For the numerical experiments presented below, I calculated the right hand side 
of p.7p using the computer algebra package MathTensor. For the integration of the 
system of ordinary differential equations (ODE's), I used a 4th order Runge-Kutta 
integration algorithm. 



The Schwarzschild-Kruskal spacetime 

We calculate initial data for the conformal Gauss gauge in the Schwarzschild-Kruskal 
spacetime using different coordinates and conformal compactifications. Resulting 
gauges are illustrated in the numerically generated conformal diagrams Fig. 13.21 and 
Fig. [331 

The physical Schwarzschild metric reads 

where fg is the Schwarzschild radial coordinate and > 2m. We transform the metric 
using retarded and advanced null coordinates 

u = t — {fs + 2™ In (fs — 2m)) , v = t + + 2m In (r^ — 2m) , (3-8) 

After the inversion r{fs) — ^ and the rescaling with 9 = r following g = 9^g we get 

g = -r^(l - 2mr) du'^ + 2dudr + da'^ , g = -r^(l - 2mr) dv'^ -2dvdr + da'^ 

We use the metric g to calculate the right hand side of (|3.7p . The coordinates u,v 
extend analytically into regions where fs < 2m. For a simulation through we use 
the retarded null coordinate u, for a simulation through the future horizon we use the 
advanced null coordinate v. 

We give initial data on the i = slice. From the transformation p.8p and the 
inversion r = we have on this slice 

r-^^l — 2mr) r^(l — 2mr) 

which leads to the induced spatial metric 

h — r dr^ H — ^ d(j^ . 

r*(l — 2mr) 

From now on, we choose m = 2 and present the calculation only for the retarded null 
coordinate u. The domain from the event horizon to spatial infinity corresponds to 
re [|,0]. We compactify h using the conformal factor (f) ~ j^^, so that h — (p'^h is 
diffeomorphic to the standard metric on the three sphere h — dx^ + sin^ x du^ , as can 
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Fig. 3.2: Schwarzschild-Kruskal spacetime in a conformal Gauss gauge. 



rs/m 



be shown by the transformation r(x) 



2(l+sinx) 

geodesic on the i = hypersurface depending on r, is given by 



The starting point of the conformal 



1 



4 In 



1 



where we suppressed the angular coordinates. The initial data for the 1-form / is given 
as /, = -^j^^^^dr. The orthogonality condition (/^^i^)* — and the normalization 



r,(l-2r, 

requirement delivers i 



\-2r. 



Initial data for the spatial frame vector reads 



K du ■ The timelike frame vector is given by eg 



(e*)i 



l-2r» 



2rl^l - 4r* 



(l-2r,)Vl-4r, 



The time coordinate at which the conformal geodesies cross is given by Tj!± = 
± — . We have ui = . J"^,^ — r. Note that lu vanishes at spatial infinity and becomes 
unbounded at the intersection of the event horizon with the initial hypersurface. 

The resulting conformal Gauss gauge can be visualized in a conformal diagram 
which depends on the choice of the spatial coordinate and the free function k. Fig. 13.21 
is the result of a numerical calculation where k = uj/{kr + 1) has been chosen in a 
neighborhood of spatial infinity such that ^+ is a straight line in the corresponding 
conformal Gauss gauge with the slope fc > 0. At the intersection of the event horizon 
with the initial hypersurface where w — > cxd, this choice would lead to k — > cxd. To 
avoid this, we choose a different n in the interior of the event horizon which reads 
K = 2r/(l — 2r). There is a smooth transition region between the domain inside the 
event horizon and the exterior region. 

Illustrated in Fig. 13.21 is the "upper right part" of the Penrose diagram for the 
Schwarzschild-Kruskal spacetime (Fig. II. 2p . The lower horizontal line corresponds 
to the hypersurface {t — 0} in the standard Schwarzschild coordinates, where also 
{r = 0}. We see that the conformal geodesies cover in a smooth way spacelike, null 
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Fig. 3.3: Schwarzschild-Kruskal spacetime in a conformal Gauss gauge. 



and timelike infinity and the domain close to the singularity. This result is robust 
in the sense that it allows variations of initial data which satisfy the orthogonality 
condition (f^i'^)* = 0, or the use of different compactifications g based on different 
coordinates. 

Another representation of the Schwarzschild-Kruskal spacetime is given in Fig. 13.31 
For this representation a spatial coordinate x has been chosen via r(x) — 2(i+sinx) ' 
The domain from the event horizon to spatial infinity corresponds to x G ['"'/S, tt]. We 
have in this coordinate 



sin^X 



2(1 + sin x)' 



dx +sin xd<J , 



2 sinx(l + sinx) 
|cosx|(2 + sinx) 



We choose k — sinx. The causal structure near spatial infinity can be clearly seen 
in Fig. 13.31 Denoted are the sets and 1+ all represented by in the one-point 

compactification. We will refer to these sets in the following sections in the context 
of the regular finite initial value problem at spatial infinity. Note that the presented 
conformal diagrams are not Penrose diagrams as the light rays are not represented 
by straight segments with 45 degrees to the horizontal. In a Penrose diagram, spatial 
infinity is necessarily represented by a point. 

We can do the calculation of the conformal Gauss gauge using other compactifica- 
tions of the Schwarzschild spacetime as well. As an example, consider compactifying 
isotropic coordinates given by the radial transformation 



p{r) = 2 (ts- m + a/ rs (r^ - 2m) j 



The Schwarzschild metric becomes 

2 



9 = 



1 - 



mp/2y 



1 + mp/2 J 



dt 



(1 + mp/2f 



{dp^+p^da"^) 



Choosing as our initial hypersurface the t — Q hypersurface and setting the spatial 
conformal factor to (f> = (? j (1 -f- mp/2f' we get h = (j)'^h = dp^ + p^da^ . 
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The initial point has — (0, p,). We set further 

The spatial frame vector on the initial hypersurface reads ei = ndp. For the choice of k 
we calculate uj = p(l + mp/2). The choice k = p leads to a similar conformal diagram 
as Fig. 13.21 The use of different coordinate systems or different compactifications for 
calculating the right hand side of (|3.7p does not affect the quality of the calculation. 

A nice feature of the conformal geodesies is that we have some global control on 
their behavior in a given spacetime. This analytic knowledge has been used to test 
the quality of numerical calculations. One test is given by the double role that the 
conformal factor plays. While f2 is some known function given by (|3.5[) in terms of 
coordinates and initial data, the set where it vanishes corresponds to the conformal 
boundary of the spacetime and has therefore a special meaning. The spatial coordinate 
value on the conformal geodesies must correspond to ^ 00 at the set {il = 0}. 

A stronger test has been made using the behavior of the frame. The 1-form dk = 
fl ff^Cf.^ is known explicitly in terms of the initial data and the coordinates via the 
relation (|3.6p . Comparing the evolution of with its a priori known form (|3.6p 
delivers a strong test for the numerical calculation. 

The Kerr spacetime 

Numerical techniques allow us to go beyond the spherically symmetric case in the study 
of the conformal Gauss gauge. As shown in this subsection and seen by the numerically 
generated conformal diagram Fig. 13.41 we can solve the conformal geodesic equations 
in the Kerr spacetime covering null infinity, timelike infinity and the Cauchy horizon. 

We take the Kerr metric in Boyer-Lindquist coordinates with m > a where m is the 
mass and a is the angular momentum of the Kerr spacetime. We make an Eddington- 
Finkelstein-like transformation using ingoing and outgoing light rays as coordinate 
lines (see [78] for the transformations). We introduce a compactifying coordinate r 
which is related to the physical coordinate f by the inversion r(f) = 1/r. For the 
outgoing case with the retarded null coordinate u we get for the Kerr metric 

/ 2m\ 2 2 -iam . 2 a j 1 

g = — 1 du Trdu dr sm v du dip 

\ rY, J r rS 

+ '^sm^ddrdLp + Y.dd'^ + ^{^^ + a^^ - Aa^ sin^i?^ sin^i? 

where 

^1 9 9 „ . 1 2™ o 

E ^ — + cos^ A== — - — + 0^. 

The metric g with respect to which we calculate the right hand side of (|3.7p is obtained 
by 5 = O^g with 9 = r as before. 

We give initial data on a i = hypersurface in the Kerr spacetime where t is the 
timelike Boyer-Lindquist coordinate. Following [43], we do the coordinate transforma- 
tion fix) = ^ .) . We rescale h with the conformal factor 

^ ' m sin x-l- V ^Ti — a- 



sin X \/ m? 



(1 — TOr)\/l + a^r'^cos'^'i) 
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Fig. 3.4: Kerr spacetime in a conformal Gauss gauge. 



The compactified spatial metric reads 

h ^ (ji^h = dx^ + sir? xd-d^ + sir? x sir? ( 1 + l + 2m/(I]r(x)) ^.^2 ^ gj^2 ^ \ ^^2 

V Ssm X 

We write the initial data using the coordinate r. We get 



mr + a?r^ cos^i? 



^(1 



^)(1 



■dr 



and 



(1 



mr)(l + a?r^ 



cos^d) 
cos^-d) 



-j2j,2 cos^i? — 2mr 



a?r^ cos -dsind 
1 + a^r^ cos'^'d 



d^. 



Note that for a — we get the formula for the Schwarzschild spacetime calculated in 
the previous subsection. 

The conformal diagram in Fig. 13.41 has been plotted with the parameters chosen 
a,s = tt/2, m — 2 and a — 1. For the "flat" ^+ depicted in the figure, the free 
function k has been chosen according to k = uj. We made this choice only to illustrate 
the resulting conformal diagram. For a numerical solution to the Cauchy problem, 
this choice is bad as the grid speed of ingoing characteristics becomes unbounded near 
J^+. It is bad also in the interior of the event horizon as the Cauchy and the event 
horizons meet on the initial hypersurface in such coordinates. It turns out, however, 
that the choice is useful in analytic studies [1401 1143] . 

In both the Schwarzschild-Kruskal and the Kerr spacetimes the conformal geodesies 
are seen to be regular in the interior of event horizons. In the Kerr spacetime they 
pass regularly through the Cauchy horizon where one needs to change the coordinate 
representation of the metric to integrate the equations further. This property of the 
conformal geodesies suggests that the conformal Gauss gauge might be suitable to 
study the inner structure of black holes and questions on stability of Cauchy horizons 
[4TI 1103] . The behavior of conformal geodesies in the interior of black holes needs to 
be analyzed in further detail before such studies can be made. As our main interest 
in this thesis lies in the treatment of the asymptotic region we do not follow these 
questions further. 
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3.2 The reduced general conformal field equations 

To use the conformal geodesies as the underlying gauge, Friedrich extended the con- 
formal field equations to admit Weyl connections. This extended system is called the 
general conformal field equations and is written for a conformally rescaled frame field 
with respect to a Weyl connection. 

The introduction of a rescaling and the transition to Weyl connections allows a more 
general conformal gauge freedom than the one studied in chapter [2l The conformal 
gauge now includes in addition to the conformal factor fi, a 1-form /. We have 
seen that in vacuum, conformal geodesies determine a conformal factor and a 1-form 
defining a Weyl connection. Friedrich used the conformal Gauss gauge based on the 
conformal geodesies to fix the extended conformal gauge freedom. He derived a reduced 
system that is equivalent to the vacuum Einstein equations and implies a symmetric 
hyperbolic system which preserves the constraints on a conformal extension. We call 
this system the reduced general conformal field equations. 

Based on the reduced general conformal field equations, Friedrich formulated a 
regular finite initial value problem near spatial infinity in which the location of spatial 
and null infinity are known a priori [61j . The construction allows one to study in 
detail properties of solutions to Einstein equations near spatial infinity (see [62l El] 
for reviews). 

The reduced general conformal field equations are written for the following vari- 
ables: a frame field e,!^, a Weyl connection T^^j, the Schouten tensor Lij = \R(ij) — 

\R[ij\ — j^Rgij, and the rescaled conformal Weyl tensor W^^^j,; = ^^jki- Note that 
the rescaling is such that if the conformal Weyl tensor satisfies the Sachs peeling behav- 
ior, the rescaled conformal Weyl tensor attains a regular limit at . We set = Cq^ 
as in 13.1.21 to construct a conformal Gauss gauge. We have as a consequence 

eo^ = 5o^ .fo = UeJ'^0, f = 0, Lqj = 0- (3.9) 
The reduced general conformal field equations read 

ViW\^, - 0. (3.10) 

The functions f2 and dk are known a priori in terms of initial data and are given by 
([33)1 and jSH). 

The numerical and analytical advantages of the system (|3.10l) are similar. We take 
a numerical point of view: 

• The system is regular for all values of the conformal factor so that no calculation 
of formally singular terms is needed that might introduce numerical instabilities. 

• The system consists mainly of ordinary differential equations (ODE's) except 
the Bianchi equation, which admits symmetric hyperbolic reductions. Besides 
enhancing the accuracy of the code and lowering the computational cost, this 
property is especially advantageous when one is dealing with complicated geome- 
tries. 
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• The location of the conformal boundary is known explicitly so that the code can 
be adapted to calculate the physical part of the conformal extension only. 

• The Weyl tensor is one of the variables. Combined with the knowledge of the 
conformal factor, this property simplifies the numerical extraction of radiation. 

• There is a wealth of analytic knowledge on the solutions of the system near spatial 
infinity which can be used to test or improve the numerical code |61l 1141] . 

Some difficulties with the system are the following 

• The equations become degenerate at the sets where null infinity meets spatial 
infinity. 

• For non-vanishing ADM-mass, the regular finite initial value problem at spatial 
infinity imposes a spherical grid topology so that simple codes based on Cartesian 
grids cannot be used. 

• The number of variables and equations is higher than in common reductions of 
Einstein equations. 

• Numerical experience with systems based on frames is sparse. 

Note that while the degeneracy of the equations at X+ is an intrinsic difficulty, the 
other items depend on the practical experience gathered in numerical calculations. In 
any case, the advantages of the system seem to counteract its difficulties. 

We will deal with the degeneracy of the equations at the set 1+ by freezing the 
evolution in the unphysical domain given by < and by choosing a suitable time 
stepping. To achieve this, a priori knowledge of the conformal factor along the confor- 
mal geodesies is very helpful. 

As to the geometry imposed by the cylinder at spatial infinity, numerical codes are 
available that can handle coordinate singularities on a spherical grid [92l I122[ 1135] . 
Indeed, we will use an infrastructure based on the Cactus framework |137j called 
GZPatchSystem that can handle different coordinate patches to cover the sphere in a 
smooth way [135(1136] . To implement a frame-based evolution system on a sphere, we 
need to deal not just with different coordinate systems but also with different frames as 
the sphere can not be covered by a single frame field. We will extend GZPatchSystem 
such that frame-based evolution systems can be solved on a spherical grid. The fact 
that the main part of the reduced general conformal field equations determining the 
geometry consists of ODE's will simplify the numerical implementation significantly. 

The high number of variables is a common feature in conformally regular field 
equations. One may think of this as the price for having a system which is regular 
for all values of the conformal factor. In the reduction based on conformal geodesies, 
the high number of variables is mainly just a memory issue as the main part of the 
resulting equations consists of ODE's that are cheap in terms of computational time. 
These issues will be discussed further in later sections. First, we rewrite the equations 
in a form suitable for numerical implementation. 
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3.2.1 Rewriting the equations 

We want to write the equations for independent evolution variables such that only 
spatial indices appear. The components of a Weyl connection F^-'^ do not have the 
symmetries of the components of a metric connection F j-*^ . The components of a metric 
connection satisfies Fy^ = —Tikj which is convenient in calculations. We therefore 
write the Weyl connection in terms of components of a metric connection and a 1-form. 
The relation of the Weyl connection components to the metric connection components 
oi g = f2^5 reads 

rA = r/, + <5/A + s^n - v^kiffi- (3.ii) 

We count 21 independent variables (/a, F^''q, F^^^). We will rewrite the equations in 
terms of these variables. By (13.91) and (|3.1ip we have 

r 

The equation for the frame components becomes by summation 

(JQ^a — JaOQ i a O^b ■ 

The equation for F^' ^ consists of three parts 
For the Schouten tensor we get 

doLab = -F/oLefa+doM^V+rfeVK'f.oa- 

All together we get 45 ODE's for the variables {ej^,fa, ^ao, ^ac^ ^aj)- The rest of the 
system is given by the Bianchi equation that plays a fundamental role. We set n — eq 
as in |68| and decompose the Bianchi equation with respect to n and its orthogonal 
component with the induced metric hij = gij + riirij = gij + r?oi'?oj- We denote by 
Qijk the totally antisymmetric tensor with £9123 — 1 and set Eijk = n^eiijk = f-oijk- 
The electric and the magnetic parts of the rescaled Weyl tensor are defined by Eij — 
hTh^Wmknin'^n' and B,j = h^h'^W^^^^^iJi'^n' where W^*^; = iW,,>„„e"\;, so we have 
EijU^ = 0, i?/ = and BijU^ = 0,-8/ = 0. Setting = hij + niUj we write the 
splitting of the rescaled Weyl tensor as 

In our case where n* — Sq we get for the electric and magnetic parts Eab — WoaOb = Eba 
and Bab = W^*o&o = -\Woacd(-''\ = Bba- Using the relation Cacdc'""^ = -25^ we get 
Woabc = ^hc^Bda- The splitting of the rescaled Weyl tensor becomes 

Wijki = -rjoiEjkm + VojEikTioi + miEjimk - VojEiiTjok - 

VOiBJ^emkl + 11f)jEi"'€mkl + ^ij"^ Bmk1]M " ^ij"^ Bmirjok — ^ijmE'"^"' Enkl ■ (3.12) 
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We use these relations to replace the Weyl tensor by its electric and magnetic parts. It 
is convenient for numerical implementation to write all evolution variables with lower 
frame indices. We get 

dofa = —fc^adoS'^'^ + LaO, 
C^orabO = —^cbO^adoS^^ + Lab — Eab, 
do^abc — { — ^dbc^aeO + Sabfd^ceO — Sacfd^beo) S'^'^ , 

— SabLcO + SacLbO — fc^abO + fb^acQ + ^ ^bcdBae^'^'^ , 
doLaO = { — ^acoLeO—dcEea)^'^'^, 

doLab = {-TacoLeb + dce,kdBaf5''^W'-doEab. (3.13) 

The rescaled Weyl tensor has 10 independent components. We use the tracefreeness 
property of its electric and magnetic parts (-Eafc, Bab) for the replacement -Ess = —En — 
E22 and -B33 = —Bii — i?22- 

The splitting of the Bianchi equation, yiW\^f, = 0, has been given in [59l[68]- We 
set Kij = hfVknj, K = h^^Kij, a' = n^Vjn\ VkEij = h^K^^K^V lEran- In our gauge 
we have 

We get the following evolution equations for the electric and the magnetic parts of the 
rescaled Weyl tensor 

CnEab + V.Bdiaebf " SE\^Kt)d + KEab - ea'%'^E,,Kdf = 0, 
CnBab-V,Ea(^aebf-^B\aKb)d + KBab-ea^''e,'^B,eKaf = 0, (3.14) 

where £„ denotes the Lie derivative along n — cq- The system (I3.14p is not yet 
ready for implementation in a code based on method of lines. The derivatives along 
the spatial frame vector fields include time derivatives because the time components 
of the spatial frame do not vanish in general, that is e^" 7^ 0. The system p. 141) 
is a homogeneous, first order, symmetric hyperbolic system for the unknown u — 
(Eab, Bab), ignoring the tracefreeness of Eab and Bab- We can write this system in the 
form 

+ Fu = 0, (3.15) 

where A^^ and F are matrix valued functions of the unknown u. For the time integration 
using method of lines we need to build Oqu — —{A'^)~'^{A°'daU + Fu). Instead of a 
direct calculation, it is more practical for the numerical implementation to calculate 
and store {A'^daU + Fu) using finite differencing and then to build —{A^)^^{A°'daU + 
Fu). This leads to less terms in the equations which reduces the time for computation 
as well as for the search of eventual errors in the code. We use the tracefreeness 
property of Eab and Bab to replace i^aa and i?33. 

The splitting of the Bianchi equation implies the following constraints 

V-E.a + ^K'^'e^^^B^^a - 
V-B,a + e:\2K\- K^^)Ebd = (3.16) 
The calculation of the constraints includes time derivatives if e" 7^ 0. 
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3.2.2 Initial data for the reduced conformal field equations 

The construction of initial data for the reduced general conformal field equations from 
a given asymptotically flat solution to the vacuum constraint equations can be found 
in [64] . We will focus on the special case of time reflection symmetric data for which 
the second fundamental form vanishes. 

Assume we are given on a three dimensional manifold S with an asymptotically 
flat end an asymptotically flat solution metric h to the vacuum Einstein constraint 
equations. The data we consider is such that (5, h) admits an analytic conformal 
compactification at spatial infinity. The analyticity of the conformal compactification 
is not essential for our studies. This assumption is made for convenience. To set initial 
data for the reduced general conformal field equations on 5 = S U {i} we need to 
calculate from h the following 55 variables: {e^°' , ej* ,rj' ^, Tj'„, fa, Lab, Lao, Eab, Bab)- 

We introduce compactifying coordinates on S and choose a conformal factor 
such that the rescaled metric h = (fPh is analytic on S and = 0, Daipli = and 
DaDp(f)\i = 2 hap, where is the Levi-Civita derivative on S defined by h. In these 
coordinates we find an adapted spatial frame to h such that h{ea,eb) = ^ab- Note 
that this frame is not unique. We can use any other frame = RaCb where Rj" is a 
rotation matrix such that the orthonormality requirement h(e'^,e'fj) — Sab is fulfilled. 
We can also use different frames on different domains and patch the solution together. 
Having chosen a frame, we calculate the spatial connection coefficients ^^^i,, the spatial 
Ricci tensor Vab and the Ricci scalar f on 5 with respect to Cq. 

For the calculation of initial data we need the following intermediate quantities 

9^ __l_ _l_ 

tab = DaDb4> ~ - Sab S" DcDd4>, Sab = rab- Ti ^ab r. 



V^bD^^flM 3 3 

(3.17) 

The derivative Da is taken with respect to Ca, so that for example Da4> — &a{<P) — 

We choose a free function k which determines the value of the time coordinate at 
and the representation of spatial infinity in conformal Gauss coordinates. Then we set 
the initial data for our 55 evolution variables on S. The frame components e" are given 
by the rescaling = k e" and we set = 0. The spatial connection coefficients Tj'^ 
are calculated from the rescaled frame Ca- In the conformal Gauss gauge, cq is chosen 
orthogonal to the initial hypersurface so the connection components TJ'q correspond 
to the second fundamental form on the initial surface. It holds r„''Q = due to time 
reflection symmetry. The 1-form / on 5 is given such that f\s = (t)~^d(j) which implies 
f\s = K~^dK. Initial data for the Schouten tensor is symmetric due to time reflection 
symmetry and is given by 

Lab = k'^ (^--^tab + ^rSab^ = Lba, LaO = 0. (3.18) 

The electric and the magnetic parts of the rescaled Weyl tensor read 

Eab = -H^^tab+^Sab) , Bab = 0- (3.19) 
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The a priori known variables for p.l3p are as given in p.Sp and (|3.6p 

^(r) = ^(i~t'^), do{r) = -2r ^, rf, = 6/9,0. 

The choice k = 1 corresponds to the point compactification at spatial infinity. 
We have seen in (|3.2p that the Weyl tensor blows up near spatial infinity with 1 /r^ if 
spatial infinity is represented as a point. As seen in (|3.19p . choosing k ^ r compensates 
for this singular behavior. This choice leads to a representation of spatial infinity as a 
cylinder depicted in Fig. 13.11 



3.3 The Cauchy problem in spherical symmetry 

A first step in numerical work with the reduced general conformal field equations is 
the spherically symmetric case. From the analytic work of Friedrich f63] we know 
that the conformal Gauss gauge covers the complete Schwarzschild-Kruskal solution. 
Our studies in 13.1.31 show that the gauge is robust enough for numerical calculations. 
Therefore, we can expect that the Cauchy problem for the reduced general conformal 
field equations in spherical symmetry with initial data from the Schwarzschild-Kruskal 
spacetime can be solved numerically without major difficulties. 

We set the frame such that ei shows in the radial direction. Due to spherical 
symmetry we can assume that the frame matrix has the form 



/ 


1 








' \ 




6,0 























v 











62^/ sin 1? j 



where the coefficients depend on {t, r) only. In spherical symmetry the reduced general 
conformal field equations simplify considerably. If we are interested only in the single 
non- vanishing component of the Weyl tensor En, we just need to solve the following 3 
coupled ordinary differential equations for r2^0' -^22 and En which decouple from the 
rest of the system 

9o^2 ^ -(r2^o)^ + ^22 + 2^-^11' 

doEii = — 3r2^o^ii- 

The degeneracy of the equations at the set 1+ is not present in this simple case. 

We calculate the initial data using isotropic coordinates as in section 13.1.31 and in 
[61] . We set r2^olT=o = and choose k = p. The initial data for the Schouten tensor 
reads 

r I _ mp 

(1 + mp/2Y 

For the electric part of the Weyl tensor we get ShIt^o = —2m. The conformal factor 
reads 

n = -, — ^(1 



(1 + mp/2)2 \^ (1 +mp/2)2 



3.3 The Cauchy problem in spherical symmetry 



66 



The one form dk becomes 

2p 



di) = ~2T—— j—j, di 



(l + mp/2)4' (l + mp/2)3' 

The solution to the Cauchy problem generates the Schwarzschild-Kruskal spacetime 
in the conformal Gauss gauge. The resulting conformal diagram is similar to the ones 
presented in 13.1.31 

To check the solution, we compare the Weyl tensor in the frame representation 
adapted to the conformal geodesies constructed in l3.1.3l in the given background of the 
Schwarzschild-Kruskal spacetime with the Weyl tensor that we get from the solution 
to the Cauchy problem described above. For this comparison, we need to calculate 
the timelike coframe <t^ ^ for the conformal geodesies satisfying the orthonormality 
relations 

e»^cr-'^='5/ and a\a\riij ^ g^.^,. (3.20) 
In a frame adapted to spherical symmetry we get 



The rescaled Weyl tensor in a conformal Gauss gauge on the Schwarzschild-Kruskal 
background is then calculated by 

On the background of the Schwarzschild-Kruskal spacetime in the conformal Gauss 
gauge using isotropic coordinates we calculate 

a<^o(-eo^eA/ + e„°(e/)^). 



p(l + mp/2) 

This background calculation agrees with the result of the numerical solution to the 
Cauchy problem up to converging numerical errors. 

The calculation carried out above is the first numerical calculation of an entire, 
asymptotically flat black hole spacetime including spacelike, null and timelike infinity 
and the region close to the singularity [154'. Difficulties related to the mathematical 
representation of spatial inflnity and numerical resolution loss in the neighborhood 
of timelike infinity prevented earlier attempts of other authors to cover the entire 
Schwarzschild-Kruskal solution (see [5T1 1123| ). 

While numerical tests of tetrad formulations in spherical symmetry can be instruc- 
tive [25] , the simplicity of the reduced general conformal field equations in spherical 
symmetry docs not allow us to draw representative conclusions for the general case. 
Still, the study shows that global numerical calculations are possible with this system. 
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3.4 Regular data at spatial infinity as a point 

We want to include gravitational radiation into our discussion in a simple setting. 
We would like to study numerically a spacetime that allows us to represent spatial 
infinity both as a point and as a cylinder. The point representation is regular only for 
vanishing ADM-mass. 

To any asymptotically flat, conformally non-flat, static vacuum solution one can 
construct conformally related initial data for a spacetime with vanishing ADM-mass 
but a non- vanishing radiation held along ^ |57J . In this section we will see how such 
data can be calculated. 



3.4.1 Radiative solutions with vanishing ADM-mass 

An asymptotically flat, static solution g to the Einstein vacuum held equations can be 
written in coordinates adapted to a hypersurface-orthogonal, timelike Killing fleld as 

~g = -v^dt^ + h, v = v{x"), h = haf3dx°'dx^, (3.21) 

where v is the norm of the Killing field that in these coordinates takes the form dt- 
The induced positive definite metric on the hypersurfaces St = {t = const.} is denoted 
by h. The Einstein vacuum field equations reduce to the static vacuum field equations 
on S = St 

Paf3[h] = ^ DaDpv, h^f^babpv = 0. (3.22) 

In addition to the asymptotic flatness condition p.ip for h we require w — > 1 as f — > oo. 

Beig and Simon show in [TS] that if the ADM-mass of an asymptotically flat, static 
spacetime does not vanish, the conformally rescaled metric h 'm a. suitable conformal 
and coordinate gauge is not just smooth, but analytic at spatial infinity. A suitable 
conformal gauge is given by 

h = (f^h with (j) = (^-—^\ ■ (3.23) 



In this conformal gauge, we have f [/i] = as a consequence of the Einstein equations 
[S7] . The rescaled metric h extends as an analytic metric to i, but the conformal data 
for the evolution equations is still singular at i. Following [57] we construct from h 
conformally related initial data hrad for a spacetime with vanishing ADM-mass but a 
non- vanishing radiation field via 

hrad^a^^h with cr=f — ^—^1 • (3.24) 

\m 1 + V J 

The conformal data constructed from (h, a) is regular at spatial infinity represented as 
a point i which implies that the ADM-mass of h^ad vanishes. If hrad is not conformally 
flat, the radiation field along null infinity does not vanish as shown in [57] and therefore, 
by the positive mass theorem, we conclude that the initial hypersurface S can not 
be complete. Spacetimes constructed from such data will probably include naked 
singularities and do not represent physically reasonable solutions. They serve, however, 
as a good testbcd to study the neighborhood of spatial infinity in the presence of 
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radiation. In this thesis we are not interested in the singularity and will confine our 
evolutions in a neighborhood of spatial infinity. 

The static vacuum field equations (|3.22p in the conformal gauge p.23p imply with 
the notation of (|3.17p and (p — where ^ = the relation [64] 

DaDbcr - - Act Jab + ct (1 - /ict) fah = tab + - fia)sab = 0. 

Using the above relation the initial data for the fields Lab and Eab p.l8|.[3rT9|) simplifies 
to 

f Kfl - 

Lab = tab, Eab = Lab- (3.25) 

CT 1 — /iCT 

3.4.2 Weyl solutions 

In this section we calculate initial data of the type described above from Weyl solutions. 
These are asymptotically flat, static, axisymmetric, vacuum solutions to the Einstein 
equations [145j . A metric from the class of Weyl solutions can be written in standard 
spherical coordinates as 

where the functions U and K arc independent, U satisfies the flat-space Laplace 

equation. At/ = = {^-^drir^ U , and K is determined from U by 

quadrature up to an additive constant [181 1128| . Due to the requirement of asymp- 
totic flatness, we have C/ — * and i^T ^ as f oo. Weyl solutions are uniquely 
parametrized by the asymptotically flat solutions U to the Laplace equation in flat 
space which can be written in the form 

oo 

C/= ^a„f~("+i)p„(cos,?), 

ri=0 

where the a„ are constants and the P„ are the Legendre polynomials. We introduce 
the compactifying coordinate r — 1/f which maps spatial infinity to the origin. The 
metric h induced by g on a {t = const.} hypersurface reads 



{dr^ + di}^ + e-'^^r^ sin^ d dip^) 



For Weyl solutions a simple choice for the conformal factor is given by 0' = r"^ ^ 
which leads to the analytic metric h' = (j)''^h = dr^ + dd^ + e~^^r^ sin^ By 

comparing with p.2ip we see that v — , so that the conformal factor (/) — ^ 

also leads to an analytic metric. To calculate initial data for the reduced general 
conformal field equations, we can use the simple metric h' or the metric h = (fP'h 
which leads to the simplification (|3.25p . The conformal relations for these metrics are 
as follows 



h' = (t>'^h, 0' = r-e 



2M-K 




(3.26) 
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The physical initial metric induced by the radiative spacetime with vanishing ADM- 
mass is given by 

-2- . A 2 f/^^ 

hrad — cr h, with cr = — tanh — 

\m 2 

We can write it also as 

hrad = a-2(00'-i)2/,' = ^'-^h', with a' = ^ = e-^r^ cosh-^ ( ^ 

(P V 2 

The advantage of using {h' , a') for calculating initial data is the simplicity of h' . The 
pair {h, a) is advantageous because in this conformal gauge we have f[h] = and the 
related simplification (|3.25p in the construction of initial data. 



The Curzon solution 

The mathematically simplest Weyl solution is the Curzon solution ,39J which is given 
by 

m w? sin^ d 

^ = --' ^ = 2f^- 

In the compactifying coordinate r = 1/f we have 

U — —mr, K ~ —-m^r"^ s\t? 

2 

For an interpretation of this solution see [18] and the references therein. 



3.5 A Cartesian implementation 

We solve numerically the reduced general conformal field equations for a massless, axi- 
symmetric, asymptotically flat, radiative spacetime with initial data from the previous 
section. The code is based on an equidistant grid in Cartesian coordinates {x, y, z} 
and does not make use of the axisymmetry of the spacetime. 

We expect a singularity in the interior. No attempt has been made to study 
questions on the nature of this singularity, such as the (non-)existence of a horizon or 
the precise blow up behavior of fields. Our interest lies in a neighborhood of spatial 
infinity including a piece of J^"*" . Our aim is to calculate the radiation field along 
and to show that it does not vanish in accordance with the theorem presented in |57| . 
The spacetime is axisymmetric and therefore we may expect that the radiation field 
in the direction of the axis will vanish. 



3.5.1 The initial data 

As we solve a frame-based system using Cartesian coordinates, we need to choose a 
suitable frame in which the initial data for a numerical implementation of the reduced 
general conformal field equations can be calculated as described in 13.2.21 For the 
presentation of an adapted frame I will use the metric h' as the formula are shorter. I 
used the frame Ca in the actual calculations involving the simplification p.25p which 
is related to by the rescaling Ca = [4>' 4''^)^'a- 
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Spatial infinity is mapped to the origin of our coordinate system. The numerical 
code based on this one-point compactification uses Cartesian coordinates which are 
regular at the origin. A frame adapted to the axisymmetry of h' can not be regular 
at the origin. To find a regular frame, we write the metric h' in Cartesian coordinates 

{x,y,z) 

h' = dr^ + + e'^^ r^d(p'^ = 



+ y X'' + + y 

We choose one of the frame vector fields to be d^- By choosing e'^ to be proportional 
to we get 



ei - y ^2+e-2/fy2 



{\-e^^)xy ^ , je^^x^+y^^ 



2 ^(a;2+y2)(,^2+g-2/fy2) " ^ + y2 V 

This frame is regular at the origin (remember that we require K\i — 0). Using cylindri- 
cal coordinates (p, z, ip) instead of Cartesian coordinates simplifies certain calculations. 
In cylindrical coordinates the metric reads h' — dp^ + dz"^ + e~^^ p^d(f^ and the frame 
above becomes 

cos (fidp , 



\/ COs2 if + e 2A' gjj^2 ^ 



•\/ C0S2 ip + 6^2^" 311^2 ip 



-K ■ a , 6'=COS(p \ , 

e sm V3 dp H , 63 = dz. 



We can do the calculations in cylindrical coordinates and then transform the data to 
Cartesian coordinates by a simple point transformation. 

The rest of the calculation is tedious but straightforward and will not be presented. 
I just write the conformal factor because of its importance for the analysis and the 
conformal boundary. For k — \ which results in the one-point compactification i ^ i'^, 
the conformal factor given by (|3.5p becomes for data from the Curzon solution 

n = a [1-^]=^ tanh^ t- n ■ 3.27 

Here, m is the ADM-niass of the static Curzon solution. Remember that the radiative 
solution that we construct has vanishing ADM-mass. 
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3.5.2 Form of J'^ and the computational domain 

The conformal boundary of the spacetime is given by the zero set of the conformal 
factor which is exphcitly known p.27p . From this formula we read the value of the 
time coordinate at depending on the space coordinates as 

_ (l~e""-)'(l + e"^'') 

We refer to the spatial variation of the value of the time coordinate along as the 
"form" of In the present compactification, the form depends on the ADM-mass 
of the static Curzon metric, m. In the code, we choose m — 1. The resulting form of 
J>^+ on the spatial coordinate domain [—1, 1] is shown in Fig. l3.5l on the a;z-plane with 
2/ = and in Fig. 13.61 on the a;y-plane with 2 = 0. We see that the form is symmetric 
on the xy-plane but not on the xz-plane. The diagrams show an awkward behaviour 
for a numerical calculation as the time coordinate of decreases while we move to 
the interior in the x and y directions. 




Fig. 3.5: J'^ plotted on the xz-plane. Fig. 3.6: plotted on the a;?/-plane. 



This behavior can be changed by a suitable choice of the free function k away from 
. As we are interested only in a neighborhood of spatial infinity, however, we will 
simply confine our analysis to a region depicted in Fig. 13.71 

Our computational domain is given by z} S [—0.2,0.2] and t e [0,0.2]. We 

restrict the analysis of the numerically generated spacetime to the domain {x, y, z\ € 
[—0.1,0.1] and t S [0,0.1]. In terms of the physical coordinate f, this domain corre- 
sponds to f e [10, oo] on the initial slice S. Note that we are also able to calculate a 
piece of which is infinitely far away from the singularity in null directions. This 
allows us to calculate the Weyl component ip4 in terms of a suitably adapted frame at 
null infinity as described below. 

In our calculations, spatial infinity is included in the numerical domain, but it is 
not on a grid point. The reason is that the numerical evaluation of initial data at i 
causes difficulties due to formally singular divisions by the conformal factor. This can 
be remedied by taking the limit analytically at i. For our purposes, however, using a 
staggered grid where i is between two grid points works sufficiently well. 
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Fig. 3.7: The domain of analysis in a neighborhood of spatial infinity. 



3.5.3 The code 

The code for the numerical solution of the reduced general conformal field equations 
with the point compactification at spatial infinity using Cartesian coordinates is based 
on the Cactus computational infrastructure [137j . I used MathTensor and Kranc 
[55] to generate the formula for the initial data and the right hand side calculations. 
For the time integration I used a 4th order Runge-Kutta algorithm. Spatial deriva- 
tives are calculated using second order accurate finite differences. For a grid function 
4>{x, y, z) — > (jjijjk, second order accurate finite differencing along a direction, say along 
the a;-direction, is given by 

dx(j){x,y,z) Dl4)i^j^k = j:^^{4>i+i,j,k- (t>t-i,3,k)- (3-29) 

The reason for using second order accurate derivatives instead of fourth order 
derivatives lies in the boundary treatment. To avoid the presumably naked singularity 
in the interior of the spacctime, we use an artificial timelike inner boundary. In our 
case, the numerical outer boundary, i.e. the outer boundary of the computational 
domain, lies in the interior domain of the spacetime. Therefore in the discussion, we 
refer to the computational outer boundary as the inner boundary. 

As the inner boundary of our evolution is not a generic aspect of our approach, 
no special numerical inner boundary treatment has been devised for our system. The 
inner boundary is present in this test case due to the naked singularity in the spacetime 
under study. In cases with regular data given on a complete hypersurface 5, we would 
not have such a timelike boundary. Therefore we are not concerned with the problem 
of boundary treatment and set the derivatives simply to zero at the inner boundary. 
This results in large numerical errors in the interior domain which propagate outwards 
faster than the grid speed of the physical characteristics because the treatment is not 
well-posed. Second order accurate finite differencing results in a slower propagation 
of errors than a fourth order accurate differencing. We will see that in the restricted 
domain of our analysis the code is second order convergent. 

In Fig. 13.71 the domain of analysis has been plotted for the medium resolution. 
The analysis has been done for two resolutions: medium and high. The medium 
resolution has 40 points, the high resolution has 80 points in each spatial direction on 
the computational domain [—0.2,0.2] which implies 20 and 40 points on the domain 
of analysis in each spatial direction. 
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3.5.4 The radiation field 



We calculate gravitational radiation represented by the absolute value of the Weyl 
tensor component ip^ at null infinity with respect to a suitably adapted Newman- 
Penrose tetrad in the x and z directions. 

In an adapted complex Newman-Penrose tetrad (Z, k, m, m) along ^+ fulfilling 
the relations (|2.26p , the complex rescaled Weyl tensor component tfj4 representing the 
radiation field is given by [32] 



(3.30) 



The advantage of knowing not only the location of J^"*" but also the conformal factor 
in explicit form simplifies the numerical radiation extraction considerably. 

The vector field I along J''^ is given by U ~ ViVl = e^^d^Q,. We can calculate 
d^Vl\j;+ a priori from (|3.27p or we can use the known field di = Vtfi + V^ri from (|3.6[) . 
I chose to use the explicit form of the conformal factor. By our gauge conditions we 
have Cq^ — and therefore also Iq is known a priori. We only need to read the values 
ej^ from the numerically generated solution to calculate the adapted frame vector field 
Z, the rest is known explicitly. 

There is an arbitrary choice involved in the calculation of ?/'4 that we fix in the 
following way. Given U = (lojhjhjs) with respect to some frame e,, we can always 
find a new frame e'^, , so that I takes the form Z^/ = (Zq' , 0, 0, Zq') with respect to the new 
frame e-, . We have Zq = Iq' — ±a/ZJ + + ^ ^^^^ ^^'^ rotation is purely 
spatial. The new frame e'^, is related to the original frame by eo = Bq and e^, — Rj'iei,. 
The rotation may be thought as consisting of two subsequent two-dimensional rotations 
resulting in 

(cosai2 — sinai2 

cosa23 sinai2 cosa23 cosai2 
sina23 sinai2 cosq!i2 sinQ!23 

Assume without loss of generality that Z2 ^ and Z3 0. The first rotation is done 
with the angle ai2 = arctan (j^^ on the 12-plane to achieve ly = and the second 




one with the angle a23 — arctan 
rotation matrix becomes 



on the 23-plane to achieve Z2' — 0. The 



h 



h/lo 



h/lo 



lo 

h/lo J 



We have used -\/Z^ + Z| + Z| = h. In this new frame, mi> can be written in accordance 
with the relations (I2.26P as m^/ — (0,0, a + ib,ia — b) where = —1 and a, 6 € R 
satisfy + 6^ = i. In our calculations, we choose a — and b ~ The rotation 
matrix is used to calculate the electric and the magnetic parts of the rescaled Weyl 
tensor in the new frame via Ea'b' = R^Rf^Eab and Ba'w = R^Rb^Bab- Using the 
splitting (|3.12p . the tracefreeness of E and B and the definition of i/'4 given in (|3.30p 
we see that ?/'4 is given in the new frame by 



V'4 = {10?{{-Evv + S2'2' - 2 B2'V) + i [BVV - B2'2' - 2 E2'l')). 



(3.31) 
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Rescaled |Psi4l2 on scri 



Radiation field on scri 





Fig. 3.8: Absolute value of V4 along ^+ in ^.^^ 3 Absolute value of ^-4 along J 

the z-direction rescaled according to reso- • ,1, j- j.- j.- j- j. 

m the x-direction representing radiat: 

lution. 



ion. 



We discuss the absolute value of ^^4 calculated from (|3.3ip along in the x and 
z directions in a neighborhood of spatial infinity on the domain illustrated by Fig. 13.71 
While we did not make explicit use of the axisymmetry of the spacetime in the code, 
we expect that it makes itself manifest in the radiation field. The radiation field in 
the direction of the axis along J'^ should vanish. This can be seen nicely in Fig. 13.81 
which shows the absolute value of 7/14 in the direction of the axis, the z-direction in our 
coordinates, along ^+ in two different resolutions. We see that the absolute value of 
'!/'4 vanishes up to numerical errors with the expected order of accuracy. The curve from 
the high resolution run has been multiplied with 2^ = 4 and lines up quite accurately 
with the curve from the medium resolution run which implies that the code is second 
order accurate in a neighborhood of spatial infinity. We also see that errors coming 
from the inner boundary destroy the convergence at a ^-coordinate distance of about 
0.08 to spatial infinity along This problem can be dealt with by devising a better 
numerical boundary treatment, but we are interested in the radiation field only in a 
neighborhood of spatial infinity and for this purpose the code is good enough. 

The radiation field along the a;-direction has been plotted in Fig. 13.91 for the high 
resolution run. To see the effect of the inner boundary on this plot, we need to go fur- 
ther in the direction of the source. We see that the radiation field grows monotonously 
as we move away from spatial infinity along J'^ . The source of radiation is presumably 
the naked singularity in the interior of the spacetime. 

We note that the simplification of radiation extraction over earlier forms of confor- 
mally regular field equations is two-fold. Firstly, we know the location of ^+ a priori 
by ()3.28p and therefore do not need to find numerically the zero set of the conformal 
factor (note though that the zero set is also known in a ^+-fixing gauge in the context 
of a hyperboloidal initial value problem [l^). Secondly, we know the explicit form of 
the conformal factor (|3.27p and not only its zero set. This allows us to write down 
analytic expressions for the derivatives dpS^\jJr^ which would have to be calculated 
numerically on J^^ if the conformal factor would be one of the evolution variables. A 
disadvantage over a ^+-fixing gauge is that the location of the extraction surface is 
not fixed on the grid. 
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3.6 Numerical implementation of frames on 5^ 

The numerical simulation presented in the previous section is based on the assump- 
tion of a regular point compactification at spatial infinity and cannot be generalized to 
include the more interesting cases with non-vanishing ADM-mass. For a generalizable 
numerical code, we need to implement the equations in a gauge in which spatial in- 
finity is represented by a cylinder. The cylinder implies a spherical geometry, because 
its construction and the resulting regularity of the conformal data distinguishes the 
direction towards the cylinder in a substantial way as will be discussed in the next 
section. The main difficulty with implementing a spherical grid geometry is that the 
sphere can not be covered by a single coordinate chart in a non-singular manner. 

One possibility to deal with this problem is to use pseudo-spectral methods which 
are promising due to their high accuracy and low memory requirements \22\ I48j . These 
techniques apply spectral methods for spatial derivatives and method of lines for time 
integration. In pseudo-spectral methods, the expected singular behavior of variables at 
coordinate singularities can be dealt with by an appropriate choice of function spaces 
which respect certain parity properties [TTl [22l |48l [88] . 

The pseudo-spectral approach does not seem to be appropriate for our problem 
for two reasons. The first reason is the nature of our initial data. We would like 
to avoid effects from the interior where presumably a naked singularity resides and 
want to study only a small neighborhood of spatial infinity. Pseudo-spectral methods, 
however, are not local methods in contrast to finite differencing methods. One might 
introduce some localization by using spectral elements instead of global basis functions. 
This roughly corresponds in the language of finite differences to high order accurate 
stencils 22, 48^. Remember, however, that in the previous section we chose a second 
order accurate finite differencing stencil to separate, as much as possible, the study 
of a neighborhood of spatial infinity from the influence of the inner boundary. It 
seems desirable for the studies we are interested in to keep the localization that finite 
difference methods offer. The second reason for not working with pseudo-spectral 
methods is that we do not only need to deal with the coordinate singularity at the 
poles but also with the singularity of the frame because can not be covered by a 
single frame in a regular way (see [17j for a treatment of S^). 

A natural solution to the problem described above is to cover the sphere in a non- 
singular manner using multiple charts on which different frames can be defined. It is 
not a coincidence that this idea which underlies the concept of a manifold has found 
its way into numerical relativity. 

To my knowledge, there are currently two approaches in numerical relativity in the 
context of an evolution problem based on finite differences that use multiple coordi- 
nate patches on a sphere. They differ by the way neighboring patches are organized. 
The " penalty method" uses touching patch boundaries over which characteristic infor- 
mation for a first order hyperbolic system is exchanged. Penalty terms drive ingoing 
modes of one patch to the outgoing modes of a neighboring patch [331 [Ml I125j . The 
"ghost-zones method" uses overlapping grids with redundant points at the boundaries 
of the patches, called ghost points, where information from the interior of a neighbor- 
ing patch is interpolated to |1161 11351 1136] . As the equations we are interested in are 
naturally written in first order symmetric hyperbolic form we can use both methods. 
In this thesis, the second approach will be followed. The reasons for this choice are 
rather circumstantial than fundamental. 
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3.6.1 Coordinates of GZPatchSystem 

The ghost-zones method has been implemented by Jonathan Thornburg in a code 
called the GZPatchSystem. The acronym GZ stands for ghost zone. I used a version of 
GZPatchSystem based on the Cactus computational infrastructure [137^ and the Carpet 
driver ^124^ . In the following we will discuss certain properties of GZPatchSystem that 
are relevant for the implementation of a frame-based evolution system. Information 
on further details can be found in the references |135|, I136j . 

From different ways to cover a sphere with non-singular coordinate charts, 
GZPatchSystem uses the "inflated cube" coordinates with 6 patches as described below. 
Another notable possibility is given by two stereographic coordinate charts [SO'. Due 
to large coordinate distortions in stereographic coordinates along patch boundaries, 
however, these seem to result in less accurate codes jll6| . 

In the following discussion, it will be convenient to think of the spheres r = const, as 
embedded into R'^ with its standard global Cartesian coordinates {x,y,z}. 
GZPatchSystem uses 3 angular coordinates denoted by /z, v and ip to cover reg- 
ularly which correspond to the rotation angles around the Cartesian coordinate axes. 
Their relation to the Cartesian coordinates is given by 

^{x,y,z) = arctan(y/z), v(x,y,z) = arctan(a;/z), Lp(x,y,z) = arctan(2//a;). 

We enumerate the patches as: 0,1,2,3,4,5 +z,+x,+y,—x,—y,—z. We use the 
pair {v, ip) on the ±x- respectively (1, 3)-patches, (/i, (p) on the ±y- respectively (2, 4)- 
patches and {i/, fj) on the ±z- respectively (0, 5)-patches such that the local coordi- 
nates on each patch are regular covering a neighborhood of the Cartesian axes. The 
internal coordinates of GZPatchSystem denoted by (p, cr) are defined on the interval 
[— 7r/4, 7r/4]. This implies the following translations for the coordinates on the patches 

P±x ^P(13) ■ {v,ip) = (pT7r/2,cr), P+,_ = P(o) : = (p,cr), 

P±v = -P(2,4) : (a*, v) = (p T 7r/2, a =F 7r/2), P^^ = P(5) : {v, fi) ^ {p - n, a - tt). 
The standard polar coordinates (p) on the unit sphere are related to (z^, /i) via 

tani?(j/, /i) ~ \J tan^ ji + tan^ tan 1^9(1/, /x) = tan/i/tani^. 

The standard metric on S'^ in polar coordinates reads ds^ — + sin ■& dpi^ . In local 
patch coordinates (p, a) this metric takes the form 

ds^ — 5 5 — - (cos^ cr dp^ — 2 sin a cos a sin p cos pdpda + cos^ p da^) . 

(1 — sin p sin a)^ 

(3.32) 

The main task of GZPatchSystem is the interpolation of a specified set of functions 
from the interior of neighboring patches into each patch's interpatch boundary ghost 
zones. We refer to this process as "synchronizing" the interpatch boundary ghost zones 
in accordance with |136j . On each patch the local coordinates are used to define the 
tensor basis. To communicate the information between different patches, one needs to 
do a a point transformation and a tensor transformation of the field variables. As we 
use a frame basis, we do not need the tensorial coordinate transformation except for 
the vector fields that constitute the frame. Instead we need to implement the rotation 
of frame-based variables. Later, we will refer to the rotation of the variables as being 
part of the synchronizing process. First we discuss the choice of frames on the unit 
sphere. 
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3.6.2 Choice of frames on S'^ 

The sphere can neither be regularly covered by a single coordinate chart nor by a single 
frame field. The coordinate singularity is independent of the singularity of the frame. 
Any frame field on the unit sphere in standard coordinates (??, ip) can be written as 

sin a . cosa 

ei=cosaa^i — -d„, e2^smad^ — -d^, 

smt* smv 

where a = a{-d, ip) is an angle of rotation. We see that at least one of the frame vector 
fields becomes singular at = 0, tt independent from a and therefore independent 
from the choice of the frame. For a = 7r/2 we recover the standard frame on the 
sphere given by = and e<^ = jj^c^i/j- In this choice, the singularity of the 
frame coincides with the singularity of the coordinates so that the frame vector field 
becomes singular at the center of the patches (0, 5). 

Given coordinates {x^,x^} and a metric g on S^, we can use the one-parameter 
rotation freedom to adapt a frame to a coordinate direction so that one of the frame 
components vanishes. Choosing arbitrarily e^^ = we get by orthonormality relations 
up to reflections 



' ^ ' ' Vdetg' 2 V.9iidetg ^ ' 

Using the above relations we can set the frame according to the metric given in local 
patch coordinates by (|3.32|) . It turns out that the standard frame {e^,e^) in the 
inflated cube coordinates is adapted so that one of the frame components vanishes on 
the patches (1,2,3,4). In terms of local patch coordinates the standard frame on 5^ is 
given by 

(1,3) 1 — COS^ V SVC? Lp (2,4) 1 — COS^ ^ COS^ 

e§ = = : d^, 

cos p sm ip 



/-, 2 • 2 / , ^3 , ^ ^ \ (2,4) 

— cos^ V sm (p — tan ip cos u Oy ^ o^^ 

\ sin V 

^ \/l - COS^ ^ C0S2 ip [ Q _i_ Q 

\ tan ip sm /i 

On the patches (0, 5) the standard frame is not adapted and becomes singular. There 
we choose a different frame given by 

(0,5) 1 .2 • 2 

ei — (1 — sm asm lyjOn, 

cos J/ 

("'■5) A ~2 1 T , • . T 

62 — — \/ 1 — sm u sm ly + sm ri tan v a, 

* \ cos ^ 

It is convenient to have the frames written in global Cartesian coordinates 

e^) = , I {xzd^+yzdy - (a;^ +y'^)d^) , , / „ i~ydx+xdy) , 

r -y/x^ + y^r yx^ + j/^ 

^1 = — / I , {xy dx + (x'^ + z^) dy - yz dz) , 62 = , „ {-z dx + xdz). 
r V a:^ + V a;^ + z^ 

(3.34) 
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Note that we have chosen only two frames on S"^, namely the standard frame (e^, e^p) 
on the (1,2,3,4) patches and an adapted frame (ei, 62) on the (0,5) patches. This means 
that we need the rotation implied by the different frames only between the 2 polar 
patches and the 4 equatorial patches in contrast to the coordinate transformations 
which need to be done between each pair of neighboring patches. The number of 
frame systems to cover S*^ should not be relevant for the accuracy of finite differencing 
in contrast to the number of coordinate charts because finite difference operators act 
along the coordinate lines and not along the integral curves of the frame vector fields. 



Rotation of the frames 

To communicate information on evolution variables between different patches, we need 

3a = Ra 



to use the rotation between the frames = R^cb- The rotation matrix can be 



written as 

j^B _ ( cos a — sin a 
y sin a cos a 

Using the orthonormality relation in (I3.20|) . the rotation matrix can be calculated by 

Rf^Rfe^^a^^^e^a^^. (3.35) 

There are different possibilities to implement this calculation and the transformation of 
frame-based grid functions between neighboring patches. We shall discuss this question 
on a simple example in the next subsection. 



3.6.3 The eigenvalue equation for the Laplace-operator on S'^ 

It is instructive to discuss the numerical implementation of frames using 
GZPatchSystem on a simple example that does not involve time evolution. We check 
the code by verifying the following eigenvalue equation for the Laplace-operator on an 
embedded sphere of radius r = ^y x'^ + + z'^ 

AY, -J1±Ay, 

We use this relation on a three dimensional domain to test the numerical implementa- 
tion of frames on a spatial slice. We calculate the Laplacian of Yim in first order and 
second order form. We write 4> := Y/m for notational simplicity. The Laplacian of 4> 
can be calculated by 

A0 = d'^' (e,"e/ + e," (a«e/)a^0 - T.^.e/a^g^) . (3.36) 

This formula includes second order derivatives of (p. Therefore we refer to the above 
calculation as the second order form. Another option is to define auxiliary variables 
4'a '■= ea(0) = e'^da4>. Then the Laplacian can be written in first order form as 

Acj, = 6-' {e.^dM - r„^,0,) . (3.37) 

The main steps in the calculation of the Laplacian of the spherical harmonics using 
GZPatchSystem in first order form (|3.37p are as follows: 
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• Fix the geometry. 

We fix the geometry by setting the frame fields (e^, e^p), (ei, 62), the radial frame 
vector field and the related connection coefficients F^'^j. It is convenient to 
use Cartesian coordinates for the input (|3.34p . The radial frame vector field is 
the same on each patch. 

Cr ^ es ^ ^{x dx + y dy + zdz). 

The connection coefficients F/^ = ^ = ^I'^i = ^2^2 = r same for 

both frame fields. The only non- vanishing connection coefficient that depends 
on the choice of our frame is 

r + ' r \/ + 

• Set the function = Yim- 

We set the function for some I and m = 0, 1, . . . , L We give the real part of the 
spherical harmonics as initial data in global Cartesian coordinates. We set, up 
to the normalization factor, 

• Transform the data to local patch coordinates. 

We let GZPatchSystem transform the grid functions to each patch's local co- 
ordinates and coordinate basis. Only a point transformation is needed for the 
spherical harmonics and the connection coefficients. For the frame fields a point 
transformation as well as a transformation of the coordinate basis needs to be 
done. These transformations are implemented in GZPatchSystem and are steered 
by the user via interface and parameter files. 

• Calculate the auxiliary variables (pa ■ 

We calculate (t>a = e^da(f> in local coordinates on each patch using fourth order 
accurate discrete derivatives as in (j2.36p . For the numerical approximation of 
da4> using fourth order accurate derivatives on each nominal grid point, the values 
of the function on two neighboring points in each direction along the coordinate 
line of x" is needed. We do not calculate the derivatives on the ghost zones. 

Fig. 13.101 illustrates a two-dimensional non-trivial stencil geometry for overlap- 
ping grids with two ghost points. The thick line is the common boundary of 
adjacent patches denoted by P and P. The solid lines are the coordinate lines of 
the patch P, the dashed lines are the coordinate lines of the patch P. Note that 
the patches share the coordinate perpendicular to their common boundary. The 
ghost zone of P is in the interior of P. For clearness the ghost zone of P has not 
been drawn. The filled small circles correspond to grid points of P, the empty 
ones correspond to grid points of P. No distinction has been made in the circles 
between nominal and ghost points. Say, we want to calculate the derivative of 
(f) in patch P on the boundary line at the grid point denoted by n in the figure. 
Beside the nominal points n — 2, n — 1, n, we need the values of ^ at n -I- 1 and 
n + 2. When initial data has been set correctly on the whole grid including the 
ghost points, the derivatives dacfi can be calculated on the nominal points. 
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Fig. 3.10: Stencil geometry with overlapping grids. 



• Synchronize the auxiliary variables. 

As we will derive the auxiliary variables for the calculation of the Laplacian in 
first order form (|3.37p . we need their correct values also on the ghost zones. These 
are not known, neither in patch P nor in patch P, but they can be calculated 
by interpolation from P. This calculation is done in GZPatchSystem using one- 
dimensional interpolation in the direction parallel to the boundary. In Fig. 13.101 
GZPatchSystem uses the empty circles in the nominal domain of P along the two 
coordinate lines parallel to the thick boundary line to calculate the values of the 
auxiliary variables on the full circles. 

• Rotate the auxiliary variables. 

For the calculation of the auxiliary variables on the ghost zone of patch P, we 
need a transformation in accordance with the rotation of the frame. We set 
ipA = Ra4>b where the rotation matrix is calculated as in (|3.35p . For this 
calculation we have two options: 

— Analytic calculation: We calculate the rotation matrix analytically using 
the explicitly known frames {e^,e^p) and (61,62) in local coordinates. 

— Numeric calculation: We use the algebraic relation (|3.35p numerically. 

The code for the analytic calculation includes switches by patch number which 
determines the coordinates and by grid side which determines the rotation ma- 
trix (each of the 6 patches has 4 sides) . The code for the numeric calculation on 
the other hand is very short. The solution of (|3.35p can be written independent 
of local coordinates, so the second option results in a code which is more efficient 
than the input of the analytic rotation matrix. This method can also be gener- 
alized to the case in which the frame components are evolution variables and the 
rotation matrix is not known a priori. We choose the numeric calculation. 

• Build the Laplacian (|3.37[) . 

• Calculate the error by 
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The calculation in low resolution is made with 20 points, the medium one with 40 
points and the high one with 80 points in both angular directions on each patch. The 
radial domain of calculation is r G [1,2]. 
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Table 3.1: Convergence factors for the calculation of the Laplacian in second order 
form. 

In Table 13.11 convergence factors are listed for the second order form of the Lapla- 
cian. They are calculated in the L2-norm for Z up to 3 on a sphere with radius r = 1.5. 
The convergence factors are independent of the sphere on which they have been cal- 
culated. For their calculation I used the code reduce by Christian Reisswig available 
from [137] including documentation. The factors have been given for pairs of patches 
that use the same coordinate system. We see a clean fourth order convergence. Note 
that for the calculation of the Laplacian in second order form, we do not need any 
synchronization. 
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Table 3.2: Convergence factors for the calculation of the Laplacian in first order form. 

Table [3?2l lists convergence factors for the first order form of the Laplacian. We see 
that the calculation is roughly third order accurate. At first sight it seems that a higher 
number of grid points results in a drop of convergence which would be problematic. 
A closer look, however, reveals that the convergence factor approaches consistently 3 
when we increase the resolution as seen for example in the case with I — 2 and m — 1. 
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The study of point-wise convergence factors using gnuplot or reduce shows that the 
convergence on the ghost points is third order which dominates in the L2-iiorm. 

The reason seems to be the interpolation of spatially derived functions and not the 
rotation. This claim can be tested by synchronizing the derivatives da(f> as covariant 
tensors instead of synchronizing the frame-based auxiliary variables (pa as functions. 
In that case no rotation is needed, but the spatial derivation da(f> is interpolated and 
derived again for the calculation of the Laplacian. One observes a similar behavior of 
the convergence factors as in Table 13.21 

In an evolution scheme, a grid function that is derived on a time step is not interpo- 
lated because synchronization takes place immediately after each evolution step in the 
MoL_PostStep schedule bin of Cactus (see [137) for a documentation on schedule bins 
of Cactus). The next subsection discusses a test case with a simple evolution system 
where we see fourth order convergence as expected. The issue with the interpolation 
of numerical derivatives will not play any role in our later discussions. 

3.6.4 Wave equation with source terms 

As a simple test for an evolution system, we solve the wave equation on a Minkowski 
background for a scalar function (j) with source / 

□0 = /. 

Writing out the equation in terms of frames we get 

The above system can be brought into first order form by using auxiliary variables 
(f'k Ek{<P)i which obey an integrability condition \7j(f>k — k't'j — 0. Take eg to be a 
timelike frame field and to be spacelike frame fields. The evolution system reads in 
terms of the frame and the connection coefficients 

eo{4>o) - 5°"^ ea{4>b) = 

For the tests we use the Minkowski background in standard spherical coordinates with 
63 taken to be the radial direction. The wave equation takes the form 

2 

9*00 = e/9i 01 -t- 61^(9201 + 9l(/'2 + 62^9202 + 9303 -f r2^l(/>l -f -03 - /, 

r 

dt(j>i = e^dicj)o + e^d2(j>i3, 

dt 02 = 62^ di 00+62^ 82 (f>a , 
9*03 = (9300, 

9*0 = 00. (3.38) 

An explicit solution to this symmetric hyperbolic system is given with the source term 
/ = -^0, by 0(t, r, 19, if) = c^it, r) Yim{'&, ^) with 

0(i,r) = -^{{r + t)e^^~^^ + (r-i)e~^^), (3.39) 
2r 

where a is some constant. The evolution consists of the following steps: 



00- 
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• Fix the geometry. 

We set the frame fields and the connection coefficients as in 13.6.31 

• Set initial data. 

The initial data consists of (j){0, x, y, z) and 0i(O, x, y, z) = {ej^d^4'{t, x, y, z))\t=o 
calculated from the explicit solution p.39p in Cartesian coordinates. Note that 
the input for the auxiliary variables 4>a depends on the frame and therefore also 
on the patch number. 

• Transform the data to each patch's local coordinate system. 

• Calculate the right hand side of the system (|3.38p and apply boundary conditions 
for integration using method of lines. 

No boundary conditions are needed in the angular directions. The boundary 
conditions we apply in the radial direction are the same as for the reduced general 
conformal field equations. We set the derivatives at the outer boundary to zero. 

• Do an evolution step. 

• Analyze the numerical solution. 

Our analysis consists of building the difference between the numerical solution 
and the analytic solution given by (|3.39p . 

• Synchronize the auxiliary variables (pi . 
This step is discussed in some detail below. 

• Calculate the right hand sides and apply boundary conditions. 

• Do an evolution step . . . 

Due to its importance in the code, we go through the steps included in the synchro- 
nization. We use the frame components to illustrate the steps involved. We write 
a two-dimensional frame field ca in angular coordinates x^ on a patch P with grid 
points n as ^e^(x„). The frame, the coordinates, the coordinate basis, and the grid 
points of a neighboring patch P will be denoted by barred indices. A synchroniza- 
tion step transforms the information from a neighbouring patch into the current patch 
^e^{xn) ^ ^e^{xn). The steps of the synchronization are (compare Fig. I3.10() 



Copying into current patch: 




- Pe 


Point transformation: 


e{x) <- 


- e{x) 


Transformation of coordinate basis: 


e A ^ 




Interpolation: 


e{xn) <- 


- e{xn) 


Frame rotation: 


CA <- 


- ^A 



The steps except the rotation are implemented by Thornburg in GZPatchSystem. 

There are three main sources of numerical errors in our evolution system that can 
be controlled to some degree: the round-off, the interpolation and the finite differenc- 
ing. The convergence factors presented below show that round-off and interpolation 
errors are negligible compared to finite differencing errors. Of course, this can only be 
expected in a high enough resolution and a higher interpolation order than the order 
of finite differencing. For the tests presented below a 6th order interpolation has been 
chosen. 
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Table 13.31 shows convergence factors for the numerical solution of the symmetric 
hyperbolic wave equation (|3.38p . The source terms have I — 2, m = 1 and I = 3,m = 2. 
The convergence factors in Table [5751 have been calculated on a sphere with r = 1.5. 
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Table 3.3: Convergence factors for the wave equation averaged over 6 patches. 

We see fourth order convergence. Tests with other choices of parameters result in 
the same qualitative convergence behavior with slightly different numerical factors. 

3.7 Implementation of the cylinder at infinity 

In the previous section we discussed the principles of a frame-based evolution code 
using overlapping grids. In this section we will apply this technique in a local study 
of spatial infinity represented as a cylinder. The spacetime under study is the same 
as in 13.51 We use a spherical grid topology instead of a Cartesian one. The code is 
designed to be usable also in studies of spacetimes with non-vanishing ADM-mass. 

We make some remarks on the structure of the reduced general conformal field 
equations with respect to the cylinder at infinity. The system consist of ordinary 
differential equations and the Bianchi equation for the rescaled Weyl tensor that implies 
a symmetric hyperbolic system as in p.l5|) . In coordinates = t, x^,A = 1,2, and 
x'^ — r we can write the Bianchi equation in the form 

{A'dt + A'-dr + A^dA)u + Fu^O. 

A rescaling with K\i = 0, dK\jr+ ^ leads to the blow-up of spatial infinity to a 
cylinder X which acts as a boundary surface to our evolution equations. This does not 
imply an initial boundary value problem as the boundary X is totally characteristic in 
the sense that A'' = on X [STJ . We get interior symmetric hyperbolic equations 
on X and no prescription of boundary data is required or allowed. The solution is 
determined uniquely by Cauchy data. 

The blow-up of spatial infinity takes care of the singular behavior of conformal 
initial data for non- vanishing ADM-mass (13. 2|) . such that we get a regular finite initial 
value problem near spatial infinity. This delicate interplay between the regularization 
of the conformal data for the field equations and the blow-up procedure of spatial 
infinity is the main reason for insisting on a spherical grid topology for our numerical 
calculations. The procedure distinguishes radial and angular directions in a substantial 
way. The smoothness of conformal data and the property that 2" is a totally character- 
istic surface is related to the vanishing of the radial frame vector field on the cylinder. 
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We can implement radial and angular frame vector fields also in Cartesian coordi- 
nates, so that the radial and angular directions are distinguished geometrically near 
the totally characteristic surface X. The points of the cylinder, however, would not 
correspond to points of our grid and this might require a more complicated numerical 
boundary treatment. The massless case studied in this section might still be numeri- 
cally stable, but the code might not be applicable in studies of physically interesting 
spacetimes. 

Another remark concerns the time evolution. As described in |62j . certain entries 
of the matrix valued function A* vanish on the set X+ where null infinity meets spatial 
infinity. This degeneracy of the evolution equations is the main difficulty in clarifying 
the open problems regarding the regular finite initial value problem near spatial infinity. 
It causes also numerical difficulties. Remember that we build in the calculation 

of the right hand side. This calculation is singular at X+ where the matrix A* is 
degenerate. The main feature that allows us to deal with this problem is that the 
location of X"*" is known a priori. A suitable choice of time stepping and freezing the 
evolution in the unphysical domain takes care of this issue, at least in the massless 
case that we study. Whether this treatment of the problem is sufficient in cases with 
non-vanishing mass remains to be seen. 



3.7.1 The initial data 

For numerical calculations with GZPatchSystem we do not map spatial infinity to r = 
as we did in the Cartesian case. Instead we map spatial infinity to a finite coordinate 
radius. The reason is that spherical coordinates are not well-defined at the origin. 
Therefore GZPatchSystem assumes that the computational domain is a shell bounded 
by two spheres with non- vanishing radial coordinate values. 

We compactify the physical Weyl solution such that the interval f e [0, oo) is 
mapped to r S (l,c») via the coordinate transformation 

'^('') = 7j = - + 1- 

r — 1 r 

It might better reflect the physical relations to map spatial infinity to some large 
radius and change the sign in the coordinate transformation so that large values of 
the compactifying radial coordinate correspond to the far field zone away from the 
source of radiation. In our case, however, we are only interested in a neighborhood of 
spatial infinity and not in the source. Therefore we have made the above choice of a 
compactifying radial coordinate. 

The conformally rescaled metric h in these compactifying coordinates takes the 
form (compare (|3.26p ) 

- ,~ /4sinh2^e^\' , , , ,K s- 2 s- 

The Curzon solution from which we calculate conformal data is given hy U ~ —m{r~\) 
and K = —^m^{r — i f sin^ d. The metric h is analytic at spatial infinity so it can be 

extended beyond r — 1. We can use the conformal factor cr = ( tanh -j)^ as before 
to generate initial data from (h,a). Another option is to use the frame adapted to 
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h' and the conformal factor a' = (r — 1)^ e^^'/ cosh(-^)- I have done the calculation 
presented below using the pair (ft,', cr'). 

A frame adapted to h' ~ {(f)' <j)^^)^'^h reads in spherical and Cartesian coordinates 

as 

ejj = ^-j- d^s = {xz + yzdy - (x^ + y^) d^) , 

r — ^ (r — 1) V 2^ + y 

^ (r-l)sm?9 ^ (r - l)Va;2 +2/2 

= 9.r = (a; 9j; + ?/ 5y + z d^)- 

The frame (e^,e;p) is used on the equatorial patches (1,2,3,4). On the polar patches 
(0, 5) around the ±z-axis, we choose a different frame (e'^, e!^) which is adapted the to 
local coordinates of the (0, 5)-patches in the sense of section [3.6.21 The frame (e'j^, e'a) 
is given in Cartesian coordinates by 



e^^x"^ +y 2 , 2 



r 



{iy^ + z'^)d^ + xydy+xzd,), 

= ((1 - e^'')xyz 9, + (e^^x^ + y^) dy + {x^ + y^)yd.), 



where / = \/y^{y'^ + z'^) + x'^{y'^ + e^^ z"^). This frame is by construction regular at 
a; = 0, y = 0. While the input of the frame into GZPatchSystem is made in Cartesian 
coordinates, the calculation of initial data using this representation is difficult because 
none of the components of the frame vector fields vanish and all coordinates {x, y, z} 
appear many times in the expressions. Instead, for the calculation of the conformal 
initial data we use the frame written in spherical polar coordinates as 

e'l — — - (cos(^sint9cost9 9^ — sin(^5„) , 

(^-1)/ 

e'n = 7 — - — — - (sin93sin?9c)^ — e^''^cos(^cos??9,^) , 
(r - 1)/ ' 



with f ~ sin cos^ ip cos^ i? + sin^ (p. Components of the frame {e'l, e'2) in this 
coordinate representation become singular at the poles given by ^? = 0, tt and these 
points are part of the patches (0, 5), but this is a coordinate singularity of the coordi- 
nates (i?, (fi) and not a singularity of the frame as discussed in section [3.6.2l on a simple 
example. We do not evaluate the above expressions on the patches (0, 5). We use it for 
the analytic calculation of the initial data as described in 13.2.21 After the calculation 
has been done, we transform the data to regular Cartesian coordinates for the input. 
Because the frame is geometrically regular on the patches (0,5), the calculated data 
will also be regular. One should reorganize the result of the calculation such that 
no formal divisions by the coordinates x and y appear (with formal division I mean 
expressions of essentially the form jj^)- Divisions by z will appear, but they cause 
no problems because on the patches (0,5) we have z 7^ 0. That such a reformulation 
can be done is due to the regularity of the data. It is a lengthy but a useful test to 
calculate the conformal data using the frames (e^,e^) and {e'i,e2) and to check that 
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the rotation of frame-based tensors agrees with the rotation of the frames on the initial 
hypersurface. 

Another point concerning the input of initial data into the evolution code is the 
divisions by (p appearing in the construction of initial data p. 181 13.19P . Though the 
data extends analytically through a straightforward evaluation of the input causes 
problems. Therefore, we evaluate the formula analytically at 

The calculation of the full data set described in 13.2.21 results in long expressions 
that will not be given as they are not relevant for our discussion. 



3.7.2 Form of ^+ and the computational domain 

Looking at the initial data for the frame e'^ we observe that the components of the 
angular frame field become singular at r = 1 which corresponds to spatial infinity 
i. The rescaling = kcJj with k ~ (r — 1) results in the blow-up of the point i to a 
sphere 1°. It leads to regular angular frame vector fields while the frame vector field 
in the radial direction vanishes on . 

A simple choice for k motivated by the regularization of conformal data would be 
K = (r — 1). Fig. 13.111 shows the resulting coordinate representation of the cylinder, 
X, and null infinity, J^+. The value of the time coordinate on is seen to depend 
on the angular coordinate in a way that seems bad for numerical calculations. It 
decreases as we move away from the cylinder in certain directions so that future null 
infinity is running backwards in grid time. This behavior leads to singularities in the 
numerical solution once the conformal geodesies approach J^+. 




S 



Fig. 3.11: The cylinder at infinity and the form of ^+ for k. — {r — 1). 

We have complete control over the form of A good choice in the asymptotic 
region that we studied also in 13.1.31 seems to be 

>^= ""^s^, ^ =fc(r-l)-f d, (3.40) 

fc (r — 1) + d 

with positive real parameters k and d. The parameter d determines the height of the 
cylinder at r = 1, the parameter k determines the slope of ^+ (see Fig. 13.12]) . 
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3.7.3 The code 

The code is based on the Cactus infrastructure |137j and uses the thorn GZPatchSystem 
|136| . We use second order finite differencing and fourth order Runge-Kutta time 
integration. No dissipation terms have been added to the evolution equations. 

We solve the reduced general conformal field equations numerically on the radial 
computational domain r e [1, 1.2]. The time domain is given by t € [0,t yr+ + At]. We 
stop the evolution after has been reached and one time step into the unphysical 
spacetime has been calculated. A reason for stopping the evolution is that the numer- 
ical solution in the unphysical region becomes singular near the set 1+ where spatial 
infinity meets null infinity. Another reason is computational efficiency. We are not 
interested in the domain beyond 

The fact that a large part of our field equations consist of ordinary differential 
equations plays a simplifying role in the implementation of the spherical grid geome- 
try. We need to implement the rotation of frame-based variables, but only for those 
for which numerical derivatives need to be calculated. These are the 10 components of 
the rescaled Weyl tensor. We do not need to implement the rotation for the remaining 
45 variables. This leads to a simplification in the synchronization step. If the evolu- 
tion of the frame components and the connection coefficients is governed by partial 
differential equations (such as for the Fricdrich-Nagy system [5^), the rotation might 
cause difficulties with GZPatchSystem because for the derivative of the transformation 
of the connection coefficients we would need to derive synchronized frame variables. 

As the geometry is not given a priori, the rotation matrix is calculated after each 
evolution step from the evolved frame components. We proceed in this calculation as 
follows: 

• After an evolution step and before the synchronization we store the components 
of the evolved frame field of the local patch in a temporary grid function. 

• We synchronize, without rotation, the 10 components of the rescaled Weyl tensor, 
and also the frame field although no numerical derivatives of the frame field need 
to be calculated. 

• After the synchronization step, the temporary grid function has the information 
on the local frame field while the synchronized frame field corresponds to the 
frame field of the neighboring patch. By using this information we calculate the 
rotation matrix on the ghost zone via p.35p . 
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• We rotate the components of the rescaled Weyl tensor as covariant two-tensors. 
After this step we can build derivatives and proceed with the evolution. 

We have two disjoint numerical boundaries in contrast to the evolution with point 
compactification at spatial infinity. The inner boundary in the physical domain (the 
outer boundary of the computational domain) is as in the Cartesian case an artifi- 
cial timelike surface. Our treatment consists of freezing the evolution at this artificial 
boundary which is not a well-posed treatment. We emphasize again that in the physi- 
cally interesting cases no artificial boundary will be present. As this inner boundary is 
not a general feature of the method we develop, we are not concerned with this issue. 

The boundary which was not present in the Cartesian case is given by the cylinder 
at infinity. The cylinder is the natural boundary to our evolution system. It is a 
characteristic surface so that no boundary data are to be prescribed. Further, as it is 
totally characteristic, no radial derivatives need to be calculated. Therefore no one- 
sided differencing needs to be implemented on the cylinder. For second order accurate 
finite differencing (|3.29p also the next point to the cylinder does not require any special 
treatment. 

3.7.4 The radiation field 

The form of in our choice of k does not depend on angular coordinates as seen in 
Fig. 13.121 Therefore, on the angular derivatives of the conformal factor vanish. 

Rescaled |Psi4l, on scri 
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Fig. 3.13: Absolute value of 7/'4 along in the r-direction representing radiation. 

Fig. 13.131 shows 1-041 calculated along ^+ in the direction of the Cartesian x-axis. 
We observe the same qualitative behavior as in Fig. 13.91 The radiation field grows as 
we move to the interior of the spacetime. The numerical values depend on the choice 
of the free function k. The parameter d that determines the height of the cylinder 
at r = 1 according to (I3.40|) is chosen to be small so that we can calculate a piece of 
before the ill-posed inner boundary treatment can destroy the solution. For the 
evolution leading to Fig. 13.131 we have chosen d = 0.02 and fc = 0.7. 

No quantitative comparison has been made between the calculation with the one- 
point compactification and with the cylinder. Certainly, more tests and experiments 
can be done in the massless axisymmetric case. It seems, however, more interest- 
ing to move on to studies of spacetimes with non-vanishing ADM-mass using further 
developed numerical methods. This will be part of future work. 
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3.8 Discussion 

The main content of this chapter has been the development of a numerical method 
that allows us to construct solutions to the Einstein equations including spatial and 
null infinity by solving a regular finite initial value problem near spatial infinity based 
on the reduced general conformal field equations. We have implemented this system 
numerically with spatial infinity represented as a cylinder. 

The method presented in this chapter is based on solving a frame-based evolution 
system with smooth inner and outer boundaries. Certain difficulties we had to deal 
with are specific to numerical calculations, such as the coordinate and the frame singu- 
larity on spheres. Other difficulties, such as the degeneracy of the equations at the set 
1+ where spatial infinity meets null infinity are shared by the numerical and analytic 
studies. It would be interesting to see whether our numerical treatment can deal with 
some mild singularity of the conformal structure at null infinity. 

We have studied a neighborhood of spatial infinity of a radiative massless space- 
time. The study of spacetimes with non-vanishing ADM-mass may require technical 
improvements of the code. We may use higher order information on the solution along 
the cylinder from analytic studies [BTJ I141j . Parallelization would allow more exten- 
sive and accurate studies of the spacetimes in question. A possibility would be to 
use parallel multi-block methods to solve evolution systems using multiple patches 

mm- 

As shown in conformal diagrams of section lSTTl the speed of outgoing characteristics 
approaches zero near spatial infinity. This property may be used to construct adapted 
finite differencing stencils. One should also consider adapting the stencil to the a 
priori known location of such that boundary grid points lie on J^~^ . For long time 
evolutions aiming at timelike infinity, one may add additional grid points dynamically 
into the computational domain such that the resolution loss due to the form of is 
compensated numerically. 

A major advantage of the conformal Gauss gauge is that the location of the con- 
formal boundary can be prescribed a priori in terms of coordinates and initial data 
by the choice of a free function k, assuming the underlying conformal geodesies are 
well-behaved. A convenient choice for numerical calculations seems to be p. 401) . Given 
a time scale and an accuracy requirement, we may choose a very large k that would 
allow an efficient and accurate calculation of gravitational radiation in the conformal 
Gauss gauge without much resolution loss in the physical part of the conformal exten- 
sion. It is an open issue which choices of k are good in the interior of asymptotically 
flat spacetimes. We have seen in numerical solutions of conformal geodesies on given 
backgrounds, that certain choices of k are bad in the interior of black holes. These 
problems seem, however, of secondary importance. The main question that one should 
study concerns the feasibility of the conformal Gauss gauge in general studies of the 
asymptotic region. If it turns out that the conformal geodesies are not well-behaved 
in the interior region, one can still devise matching methods or choose another gauge 
in a domain where the conformal factor has been set initially to unity. One can only 
hope to answer the question on the feasibility of the conformal Gauss gauge gradually. 
Possible steps for further research are outlined in section l4?2l 



Chapter 4 

Summary and Outlook 



In this thesis I took a conformal approach to the numerical calculation of asymptoti- 
cally flat spacetimes. To avoid certain conceptual and technical deficiencies of current 
numerical methods discussed in the introduction, I devised and studied new numerical 
methods which have, as expected, difficulties of their own. 

In the summary I discuss the achievements of the thesis and point out some open 
problems of the methods introduced. In the outlook I try to give a glimpse of the 
possibilities opened up by the methods assuming their problems can be solved. The 
thesis ends with remarks on the idea of conformal infinity. 

4.1 Summary 

The question that has been raised in the preface concerns the numerical methods for 
calculating asymptotically flat spacetimes, not the statements that can be made on 
such spacetimes. Accordingly the main content of this thesis does not lie in numerical 
studies of the solution space to the Einstein equations but in the development of 
new methods for future studies of the solution space. I tried to motivate the search 
for new methods in the introduction by stating some problems of currently available 
numerical methods. The results of this thesis have been presented in two chapters, 
each concentrating on a different asymptotic domain: null infinity and spatial infinity. 

Chapter 2. Null Infinity 

• Construction of a ^-fixing gauge in spherical symmetry, 12.21 

This study refutes claims made at various places (for example in [29]) that the 
conformal approach necessarily leads to a resolution loss in the physical space- 
time. The explicit examples written in a ^-fixing gauge in spherical symmetry 
show clearly that the allegedly necessary resolution loss is a property of a bad 
coordinate choice and not of the conformal approach itself, a fact sufficiently 
supported by existence results. The examples also serve as a testbed for new 
ideas in the conformal approach. 
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• The development of a numerical method that includes null infinity in the com- 
putational domain to solve a hyperboloidal initial value problem for the Einstein 
equations in a general wave gauge, 12.31 12.41 

A basic difficulty with the equations used in the conformal approach has been 
that they are very different from common formulations of the Einstein equations. 
Experience gained within the standard approach could so far not be applied in 
the conformal approach, and the clean framework that the idea of conformal in- 
finity yields could not be used in conventional numerical relativistic calculations. 
Motivated by this difficulty, 1 constructed a method to include null infinity in 
the computational domain using a common reduction of the Einstein equations. 

Assuming hyperboloidal initial data whose maximal development admit a smooth 
conformal boundary, I showed that a J^-fixing gauge can be constructed during 
the solution of a hyperboloidal initial value problem using a suitable coupling of 
the conformal and the coordinate gauge freedom in a general wave gauge. This 
result opens up the possibility to attach smoothly a conformally compactified 
asymptotic region, where conformal techniques are applied, to successful numer- 
ical calculations of the interior domain in a direct way, using the same set of 
variables and the same set of equations. 

• A numerical test of the introduced method in spherical symmetry, 12.51 

The main difficulty with the aforementioned method is the appearance of for- 
mally singular terms in the evolution equations. Within the class of spacetimes 
that admit a smooth conformal boundary and in a ^-fixing gauge, each of these 
terms attains at null infinity a regular limit. This limit can, in principle, be cal- 
culated by numerical methods. One can not assert without extensive numerical 
studies, however, that such a calculation will, in practice, result in a numeri- 
cally stable code for highly dynamical spacetimes. As a first test I studied the 
spherically symmetric case in a stationary gauge using simple numerical meth- 
ods. One could observe, in this special case, that a numerical treatment of the 
outer boundary based on extrapolation of the variables allows the calculation of 
a piece of null infinity. 

Chapter 3. Spatial Infinity 

The hyperboloidal approach does not allow us access to spacetimes in their entirety. 
To make statements on global properties of spacetimes one would like to get access to 
the conformal boundary including spatial infinity. The main problem in this context 
is the singular behaviour of conformal initial data at spatial infinity for non- vanishing 
ADM-mass. This problem was treated by Friedrich who formulated a regular finite 
initial value problem near spatial infinity using the reduced general conformal field 
equations. The construction allows a smooth extension of conformal data with non- 
vanishing ADM-mass through spatial infinity which is blown up to a sphere on the 
initial hypersurface. The development of the data results in a representation of spatial 
infinity as a cylinder. Currently the only system that allows a numerical study of 
spatial infinity including a piece of null infinity is given by the reduced general confor- 
mal field equations. The following steps have been taken in the course of a numerical 
implementation of this system. 
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• Numerical construction of a conformal Gauss gauge on the Schwarzschild-Kruskal 
and the Kerr spacetimes, 13.11 

The reduced general conformal field equations are based on the conformal Gauss 
gauge. Though it has been shown analytically that the conformal Gauss gauge 
is well-behaved in strong field regions in spherical symmetry [63j . the general 
behavior of the underlying conformal geodesies is unknown. Besides, it can not 
be guaranteed that analytical methods work well in numerical calculations. To 
test the numerical feasibility of the conformal Gauss gauge I reproduced the 
analytic construction of a conformal Gauss gauge in the Schwarzschild-Kruskal 
spacetime covering the entire spacetime in a smooth way. Going beyond avail- 
able analytical studies 1 found out numerically that one can also cover the Kerr 
spacetime using conformal geodesies including null infinity, timelike infinity and 
the Cauchy horizon. 

• Numerical simulation of an entire, asymptotically flat, black hole spacetime. 13.31 

As a first test of the numerical feasibility of the reduced general conformal field 
equations, I solved a Cauchy problem based on this system in spherical sym- 
metry with data from the Schwarzschild-Kruskal spacetime. The equations in 
spherical symmetry become a system of ordinary differential equations. This is 
a major simplification so that no statements on the applicability of the system 
in the general case can be deduced from this study. It should still be noted that 
this is the first numerical calculation of an entire, asymptotically flat black hole 
spacetime including spacelike, null and timelike infinity and the region close to 
the singularity. 

• Numerical calculation of a radiative spacetime including spatial infinity and a 
piece of null infinity, 13.51 

In order to study cases with non-vanishing gravitational radiation, I wrote an 
evolution code for the reduced general conformal field equations based on a three 
dimensional Cartesian grid using the software infrastructures Cactus |137j and 
Kranc f83j. This code can only deal with the one point compactification at 
spatial infinity which results in singular conformal data for non- vanishing ADM- 
mass. Therefore I calculated asymptotically flat, axisymmetric initial data based 
on studies by Friedrich which have vanishing ADM-mass but whose development 
has a non- vanishing radiation held. I analyzed the numerical development of this 
data and showed that the radiation content of the spacetime does not vanish but 
grows towards the interior where presumably a naked singularity resides. 

• The development of a three dimensional code to solve Cauchy problems for frame- 
based evolution systems with smooth inner and outer boundaries, 13.61 

I wrote a numerical code within the Cactus infrastructure |137j using a Cactus 
thorn written by Thornburg ^136j which implements spherical coordinates using 
multiple patches. I extended this infrastructure to include frame rotations so that 
Cauchy problems for frame-based evolution systems can be solved numerically. 

• Numerical studies of spatial infinity represented as a cylinder, 13.71 

In the representation of spatial infinity as a cylinder, the equations degenerate at 
the set where the cylinder meets null infinity. This difficulty has been dealt with 
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by freezing the evolution in the unphysical domain and choosing a time stepping 
that avoided an evaluation of the equations at the critical set. An essential 
property of the conformal Gauss gauge which allows us such a treatment of the 
problem is the a priori knowledge of the conformal factor in terms of initial data 
and grid coordinates. I suggested a convenient choice for the remaining gauge 
freedom in the conformal Gauss gauge for numerical calculations. I did numerical 
studies of the cylinder in the massless, axisymmetric case which is special in the 
sense that the spacetime under study allows the one-point compactification. 

The main open problem with this approach is related to the degeneracy of the 
equations at the critical set. Extensive studies and sophisticated numerical tech- 
niques are required before this approach can be applied to study interesting 
mathematical and physical questions concerning the asymptotic behavior of so- 
lutions to Einstein equations. 

The studies show that the a priori knowledge of the conformal factor is a very 
convenient feature in numerical studies. It seems essential that the conformal boundary 
and also the conformal factor should be controlled to some degree in a numerical 
calculation. One can expect that this feature will play an important role in the future 
development of the conformal approach. 

4.2 Outlook 

The discussed methods suggest a wide range of prospects to study the solution space 
to the Einstein equations in a way that has not been possible before. Some suggestions 
are listed below. 

• A simple but illustrative experiment would be to solve a wave equation on the 
Schwarzschild and the Kerr spacetimes written in a ,^-fixing gauge using a three 
dimensional code that can handle spherical grid topology. This might allow an 
accurate study of quasi-normal modes and tail behaviour. 

• The method of solving an hyperboloidal initial value problem for the conformally 
transformed Einstein equations in a general wave gauge by prescribing the coor- 
dinate representation of a conformal factor should be tested in three dimensions 
for dynamical spacetimes. The basic infrastructure required for this test is al- 
ready available ^106] so that one can concentrate on the main intrinsic problem 
of the method, namely the numerical calculation of the formally singular terms 
in the course of the evolution. 

• The idea of using compactifying coordinates without transforming the metric 
seems to work on Cauchy hypersurfaces [113) . although this is very awkward to 
study radiation as discussed in ll.4.21 It would be interesting to try the coordinate 
compactification technique using a hyperboloidal foliation. 

• One should construct hyperboloidal initial data for a radiative spacetime ad- 
mitting a smooth conformal boundary. There are many interesting and open 
questions related to the interpretation of hyperboloidal data that may be stud- 
ied by evolving them. 
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• A straightforward test of the numerical techniques presented in chapter 3 would 
be to evolve conformal data corresponding to static or stationary spacetimes with 
non- vanishing ADM- mass in a neighbourhood of spatial infinity and to check that 
their radiation field vanishes. 

• One should construct solutions with a complete initial hypersurface so that no 
artificial inner boundary needs to be introduced. The question how to construct 
numerically initial data for spacetimes which admit a smooth conformal bound- 
ary is open. Data whose development leads to some mildly singular behavior at 
null infinity can be constructed by solving numerically a Lichnerowicz equation 
for the initial conformal factor [43j . To numerically calculate conformal initial 
data for the reduced general conformal field equations one needs to build numer- 
ical divisions of the data by the conformal factor 13.2.21 This might not to be 
an easy task, but it seems worth the effort. A numerical study on smoothness 
properties of solutions near null infinity has not yet been done. A strong inter- 
action between mathematical and numerical studies in this question might give 
new impulses in both directions. 

• Another interesting study related to the asymptotic behavior of solutions to 
the Einstein equations would be to evolve data which might results in different 
degrees of smoothness at future and past null infinity as suggested in |142j . 

• The setting with the cylinder at spatial infinity would allow us to construct 
solutions in which pure gravitational radiation collapses to form a black hole as 
suggested in [121 [l3] (see also the discussion in 65J). 

• One might try to use use the conformal Gauss gauge in a metric formulation 
Einstein equations. This might lead to a simplification in the equations and 
require less variables. 

• The implementation of a frame-based evolution system should allow us to im- 
plement the Friedrich-Nagy system [68]. An interesting study with a successful 
implementation of the Friedrich-Nagy system and the reduced general conformal 
field equations would be to evolve a simple but radiative spacetime both with 
and without an artificial, timelike outer boundary. This would allow us to com- 
pare the approximation given by the timelike outer boundary at a finite distance 
away from the source with the idealization given by conformal infinity such that 
systematic errors can be estimated for different kinds of boundary data. Studies 
along these lines have been made without numerical access to null infinity using 
solutions for which accurate analytic expectations can be calculated [1061 1118j . 

One of the motivations behind the work of this thesis has been to develop techniques 
to calculate numerically entire, asymptotically flat, radiative spacetimes. This still 
remains an important goal for future work. A combination of numerical calculations 
near spatial infinity and null infinity might allow us to accomplish this task. An 
achievement of this goal would deliver a starting point for extensive studies of the 
solution space to the Einstein equations using numerical methods based on rigorous 
analysis. 
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4.3 Concluding remarks 

A general critique concerning the concept of null infinity and its role in physical ques- 
tions holds that null infinity is "too far away". After all, the Earth where we are doing 
our measurements is not infinitely far away from the sources and we do not move along 
null geodesies in contrast to an observer along null infinity. 

In the Schwarzschild spacetime, we have a length scale m at our disposal, so we 
can discuss astrophysical distances in numerical calculations. To give a notion for 
the scales in the conformal diagram Fig. 13.21 three curves with constant Schwarzschild 
radius = 4m, 12m, 300m have been plotted. The curve at Schwarzschild radius 300m 
corresponds to about 1800 km for a black hole with four solar masses but it can hardly 
be distinguished from ^+ in the conformal diagram. Even if we consider super-massive 
black holes, the distance corresponding to 300m is incomparably small contrasted to 
the thousands or millions of light years that separates us from the astrophysical sources. 
In this representation, the numerical effort for simulating the region from 300m to J^^ 
is very small whereas in the standard approach putting the outer boundary from 300m 
to 1000m costs considerable effort in terms of numerical techniques and computational 
sources. We note that this behavior depends on the conformal gauge. One has the 
freedom to set $7 = 1 initially on a given domain if more resolution is desired in that 
domain. Further, we have seen in chapter 2 that a positive smooth function for the 
coordinate representation of the conformal factor may be prescribed in the context of 
the hyperboloidal initial value problem such that certain spacetime domains can be 
emphasized for the calculation while still compactifying the asymptotic region. We 
may conclude that the conformal compactification technique is not only promising for 
an unambiguous outer boundary treatment and radiation extraction, but also for a 
computationally efficient numerical code to deal with the asymptotic region. 

If some astrophysically motivated length scale can be determined, as suggested for 
example in [37], such that an observer "nearly at rest in the frame of the isotropic 
cosmic microwave background radiation" at that distance from the source can regard 
the source as isolated, this length scale will naturally be given in astronomical units. 
The example above shows that a timelike curve at such a distance away from the 
source is almost indistinguishable from in the conformally compactified picture 
while the standard approach does not even allow the discussion of such length scales. 
Therefore, the arguments presented in [37j can as well be considered as supporting the 
notion of conformal infinity. After having calculated the entire spacetime, we can still 
discuss the calculation of a field which represents radiation (having a proper limit with 
respect to an adapted tetrad at null infinity to which we would have direct access) on a 
timelike surface arbitrarily far away from the source of radiation. Within the standard 
approach, however, it seems hardly possible to calculate an astronomical domain for 
an isolated system. 

In this context, I would like to emphasize that the idea of conformal infinity is an 
idealization, in contrast to an approximation. The difference can be illustrated on the 
idea of a number, say a real number. We do not have physical access to real numbers. 
Results of measurements are in a sense fuzzy. They are given in terms of rational 
numbers with error bars regarded as approximations to real numbers. Physics without 
real numbers would be, however, very inconvenient. Most physicists would agree that 
the notion of a real number is a useful one. The idealization of a real number allows a 
clean modeling in various fields of science. Similarly, the concept of conformal infinity 
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allows a clean modeling in general relativity. 

In |37j it is stated that the conformal approach implies "introducing an arbitrary 
amount of empty space not reflecting astronomical reality" . A similar statement can 
be made concerning real numbers by claiming that they imply an arbitrary amount 
of accuracy beyond the Planck scale not reflecting physical reality. The question 
whether real numbers are real in an ontological sense belongs, however, to the realm 
of philosophy. If an idealization allows a clean modeling of the 'actual reality' we can 
observe and if it is of practical use, than it is a viable physical viewpoint to use the 
concept while being aware of the limits of its application. 

Another argument claims that there is not much motivation for using the idea 
of conformal infinity, especially regarding the impressive milestones achieved recently 
within the standard approach. While a case can be made that the deficiencies of 
the standard approach are not relevant for the detection of gravitational waves as the 
accuracy of the waveforms seems to be enough for detection, there is also the viewpoint 
holding that interesting physics will be hidden behind the first couple of digits in the 
observational data. Having predictions which are accurate only up to a couple of 
percents due to systematic errors will most probably not be enough to answer many 
interesting questions in the prospective field of gravitational wave astronomy. One 
should also not forget that beside the effort of calculating waveforms there remains 
the task of exploring the solution space. Here it is unclear how useful the conventional 
approaches are if global questions are concerned. 

If an objection is to be made against the conformal approach, I believe, it should be 
made not on conceptual but rather on practical grounds. Whatever level of conceptual 
clarity, geometric appeal or computational efficiency the conformal approach promises, 
in the end, its application in numerical relativity must be judged after its practicability 
in numerical studies of the solution space to the Einstein equations. 

It seems to me that the most important task for future work in the context of 
our discussion is to devise from the conformal approach a practical tool for numerical 
calculations of asymptotically flat spacetimes. This could not be achieved yet for the 
calculation of highly dynamical, asymptotically flat spacetimes. Unfortunately, the 
results presented in this thesis also do not allow the conclusion that the suggested 
methods can be used in general calculations. I believe that this task can only be man- 
aged if insights that have been gained in the standard approach can be conveniently 
combined with the calculation of the asymptotic region. The ideas presented in this 
thesis, especially the a priori knowledge of the conformal factor in numerical calcula- 
tions, seem promising to accomplish this task. One may hope that future work will 
allow us to construct from the beautiful idea of conformal infinity a practical tool for 
numerical calculations. 



Appendix A 

Calculation of causal diagrams 



In this appendix, we present how Penrose diagrams hke Fig. ll.2l or Fig. 12.71 and radial 
Ught rays on the grid hke in Fig. 11.51 or in Fig. 12.51 for a given Schwarzschild-Kruskal 
background are calculated. We omit angular dimensions. 



A.l Penrose diagrams 



The two-dimensional Schwarzschild metric in standard coordinates is given by 

1 - — dt^+il - — \ dr?. 



The coordinate singularity at the event horizon fg = 2m can be removed by introducing 
advanced and retarded null coordinates 

V = i + r^, u = i — where = Vg + 2m In [ — 1 1 . (A.l) 



2m 



The Schwarzschild metric in these coordinates takes the forms 



/ 2m \ / 2m \ 

gs = - [l - — j dv'^ + 2dvdrs, c/s = - [l - — j dii'^ ~ 2dudrs- (A. 2) 

These line elements can be analytically extended beyond the event horizon to r > 0. 

We introduce the Kruskal-Szekeres coordinates for > 2m by 1^ = ^v/4m ^ q ^^^-^^ 
U ~ _g-«/4m ^ rpj^g nietric written in double-null coordinates becomes 



5, = - f 1 - — ^ dudv = e-'=-/2™ dUdV, 

V rs J rs 



where = fs({t, v) — fs{U, V) respectively. The Kruskal-Szekeres coordinates can be 

analytically extended to [/ > und < 0. The Schwarzschild coordinates are given 
by the implicit relations 

\2m r U 
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Compactifying coordinates can be introduced by U ^ arctan U and V ~ arctan V. 
These are defined in the range (— 7r/2 < U,V < tt/2). We define new time and space 
coordinates by setting 

t^liV + U), r^^{V-U). (A.3) 

In the Penrose diagrams, the i-axis runs vertically and the r-axis runs horizontally. 
Summarizing the coordinate relations, we have tanT^ — e^/^™, tanC/ = ^e^"/**™, 
and 

tan Uta.nV = -e^^/^" ( ^ l] , = ^e'/^™. (A.4) 

The compactifying coordinates {V, U) have been used to plot the Penrose diagrams of 
the Schwarzschild-Kruskal spacetime in Fig. 11.21 or the extended Schwarzschild space- 
time in Fig. 11.41 



A. 1.1 Ingoing Eddington-Finkelstein coordinates 

We introduce the coordinate f — v ~ fs- The metric becomes 

/ 2ni\ , 9 Am , , / 2m\ , n , . 

g, = -il- —j dr^ + —dfdr, ^ \ ~r J -* 

These coordinates are sometimes referred to as ingoing Eddington-Finkelstein coordi- 
nates in the numerical literature. The spatial surfaces that are used in a simulation 
based on ingoing Eddington-Finkelstein coordinates are given by the level sets of the 
function 

f (f, Tg) = V — = i + r^, — = i + 2m In ( — 1 

\2m 

To plot the surfaces f = const, in the Penrose diagram of the extended Schwarzschild 
spacetime in Fig. 11.41 we use (|A.4p and write (see also [98J) 

V{rg) = arctan (e^/^'^e^^/^™) , C/(f«) = arctan J^-g-^/^^e^^/^" ~ ) ' 

The Schwarzschild coordinate fg is regarded as a parameter along the curves 
{V{fs),U{'rs)) with constant values of f. We transform the above expressions using 
the compactifying time and space coordinates {t,r) from (jA.Sp . The resulting curves 
are plotted on the (t,r)-plane with Mathematica using the function ParametricPlot. 

For plotting the hypersurfaces fg = const, we can write a similar parametric rep- 
resentation. Using the relations (jA.ip and (|A.4p . we write 



UnV = e'/^^e^^/^^W — - 1, tanf/ = -e-*/4mgf,/4m / _ ^ 
V 2m V 2m 

We regard i to be the parameter along the curves {V{i), U{i)) with constant values of 
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A. 1.2 Constant mean curvature foliation 



In Fig. ll.lOl and Fig. 12.71 we plotted constant mean curvature surfaces in the extended 
Schwarzschild spacetime using compactiiying coordinates. The properties of these 
surfaces have been discussed in f96| . 

Surfaces in the Schwarzschild spacetime that are dragged along the timelike Killing 
vector field df can be written as level sets of a function ^{t,fs) ^ t — h{fs). Here, 
h{fs) is referred to as the height function. To find a CMC-slicing of the Schwarzschild 
spacetime one requires that the mean extrinsic curvature K of the surfaces $ = const, 
is constant. We get 



i:2 



For K — const, this can be integrated once with an integration constant C. A following 
algebraic manipulation results in the differential equation 



h' = 



~3~ 



C 



(A.6) 



where P{fs) is as defined in 12.2.21 We can not integrate this equation explicitly 
to get h(fs) in closed form but we can integrate it numerically. A difficulty is the 
numerical integration for the function h at the horizon for fs — > 2m. The integral can 
be calculated in the Cauchy principal value sense [96]. We build the integral in the 
regions > 2m and fg < 2m respectively via 



Kfs) 



h'{x) dx, 



and 



2m+e 



h'{x) dx, 



where e is a positive small constant. We subsequently match the resulting curves 
together at the event horizon. The embedding of the foliation is calculated as follows. 
We define the embedded surfaces by level sets of a function which can be written as 

<^ — t — h{fs) — t + — — h{fs) = V — r^t — h{fs). 

so that V = $ + (f» + h{fs))- We use the relations (|A.4p and write 



tany 

tan U 



g<I>/4mgrs/4m 



2m 



1 e 



h(rs)/4m 



_g — <['/4mgrs/4m 



2m 



1 e 



-h(rs)/4:m 



The numerical integration gives us a list of pairs of (fg, h^fs)). The surfaces corre- 
sponding to different values of $ can be plotted with Mathematica using the function 
ListPlot. 
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A. 2 Null rays on the grid 

To visualize the causal structure on a numerical grid for a spherically symmetric space- 
time in a stationary gauge we calculate the ingoing and outgoing radial null rays with 
respect to grid coordinates {t,r). The metric can be written as 

g — gttir) dt^ + 2gtr{r) dtdr + grrir) dr^ . 

A parametrized null curve t(r) fulfills the relation 

fdtV ^ dt 

For gtt ^ 0, we can calculate the curve t{r) by integrating 



t{r) = / ^"-v^tr ^^f^ (A 7) 

J gtt 

The positive sign corresponds to ingoing null rays. In the following we discuss the 
visualisation of the causal structure on the grid for the Schwarzschild spacetime based 
on different coordinates. 

A. 2.1 Ingoing Eddington-Finkelstein coordinates 

For Fig. ll.5l we use the form of the metric in ingoing Eddington-Finkelstein coordinates 
(|A.5P and integrate (|A.7|) to get 

Tin{fs) = -fs + Cm, Tout{f s) = ?s + Am\n\f g - 2m| + Cout, 

for ingoing and outgoing null rays. We plot different light rays by varying the constants 
Cin and Cout To include the excision region in Fig. 11.51 we plot the null rays inside and 
outside the event horizon separately and match them together. 

A. 2. 2 Coordinate compactification at spatial infinity 

To generate the plots in Fig. 11.61 and Fig. 11.71 we introduce a compactifying radial 
coordinate in (jA.Sp by 

r 1 

rs{r) = -j , dfs = Y^ ^ 

1 — r (1 — ry 

The domain G [3/2, oo) corresponds to r €E [3/5,1] and the event horizon is at 
r = 2/3. The metric (|A.5[) becomes 

2m, A , , 4m , , / 2m, A 1 , , 

1 (1 -r) rfr^ + — -dTdr+ 1 + (1 - r) --r dr^ . 

r J ryl — r) \ r / (1 — rj^ 

This metric is singular at r = 1 as expected. Integration of (jA.7p results in long 
expressions that can be written in explicit form. Note that the ingoing null rays in 
Fig. 11.71 that seem to come in from spatial infinity start at a small vicinity of spatial 
infinity. 
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A. 2. 3 Cauchy-Characteristic Matching 

The left part of Fig. 11.91 has been plotted using the timelike ingoing Eddington- 
Finkelstein coordinates as described above. For the right part we write the Schwarz- 
schild metric in the advanced null coordinate u as in (|A.2[1 and subsequently introduce 
a compactifying radial coordinate r via = r/(l — ?■) to get 

2m, A , ? 2 



Qs = — I 1 (1 — r) ] du — dudr. 

\ r J 1 — r 

The light cones in the parametrization u(r) are then given by solutions to the differ- 
ential equation 

2m , ,\ / du\^ 2 dii 



r ^ ^ J \dr ) ^ 1 — r dr 
The solutions are given by 

~ / \ nil I 4mln|r — 2m(l — r)| , 

Uin{r) = 2 In 1 - r — — h c™ Uout[r) = Cout- 

1 + Im 

Note that we need these light rays only outside the event horizon in the matching 
region. 



A. 2. 4 Constant mean curvature slicing 

Fig. 11.111 and Fig. 12.51 have been plotted using the form of the Schwarzschild metric 
given in (|2.13p . For this metric we can not integrate (jA.7[) to get a closed form for the 
ingoing and outgoing light rays. One needs to calculate the integrals numerically. 
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